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ABSTRACT 

Using cosmological simulations, we address the properties of high-redshift star¬ 
forming galaxies (SFGs) across their main sequence (MS) in the plane of star-formation 
rate (SFR) versus stellar mass. We relate them to the evolution of galaxies through 
phases of gas compaction, depletion, possible replenishment, and eventual quenching. 
We find that the high-SFR galaxies in the upper envelope of the MS are compact, 
with high gas fractions and short depletion times (“blue nuggets”), while the lower- 
SFR galaxies in the lower envelope have lower central gas densities, lower gas fractions 
and longer depletion times, consistent with observed gradients across the MS. Stellar- 
structure gradients are negligible. The SFGs oscillate about the MS ridge on timescales 
^ 0.4 tHubbie ("^ 1 Gyr at z ~ 3). The propagation upwards is due to gas compaction, 
triggered, e.g., by mergers, counter-rotating streams, and/or violent disc instabilities. 
The downturn at the upper envelope is due to central gas depletion by peak star 
formation and outflows while inflow from the shrunken gas disc is suppressed. An 
upturn at the lower envelope can occur once the extended disc has been replenished by 
fresh gas and a new compaction can be triggered, namely as long as the replenishment 
time is shorter than the depletion time. The mechanisms of gas compaction, depletion 
and replenishment confine the SFGs to the narrow (±0.3 dex) MS. Full quenching 
occurs in massive haloes (Mvir > -^o) and/or at low redshifts {z < 3), where 

the replenishment time is long compared to the depletion time, explaining the observed 
bending down of the MS at the massive end. 

Key words: cosmology — galaxies: evolution — galaxies: formation — galaxies: 
fundamental parameters — galaxies: quenching 


1 INTRODUCTION 

Observations of the galaxy population spanning the last 
12.5 billion years of cosmic time have revealed a picture 
in which the majority of star-forming galaxies (SFGs) fol¬ 
low a relatively tight, almost linear relation between star- 
formation rate (SFR) and stellar mass (M*), also known as 
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the “main sequence” (MS) of SFGs (e.g., Brinchmann et al. 
2004; Noeske et al. 2007a, b; Daddi et al. 2007; Elbaz et al. 
2007; Salim et al. 2007; Whitaker et al. 2012; Speagle et al. 
2014; Pannella et al. 2015). The SFR increases with M* as a 
power law (SFR oc with a ~ 1) over at least two orders 
of magnitude (~ 10® — 10^^ Mq). Several studies have found 
that the SFR towards the highest masses (M* > 10^^ Mq) 
falls systematically below the value expected for a simple 
power law relation, effectively lowering the high mass slope 


(c) 2016 The Authors 


2 S. Tacchella, A. Dekel, C. M. Carollo, et al. 


of the SFR—M* relation towards lower redshifts (Rodighiero 
et al. 2010; Elbaz et al. 2011; Whitaker et al. 2012; Magnelli 
et al. 2014; Whitaker et al. 2014; Schreiber et al. 2015). The 
most noticeable feature is that the MS relation at any given 
redshift shows a rather small scatter of (tms ~ 0.2 — 0.3 dex 
(Noeske et al. 2007a; Whitaker et al. 2012; Speagle et al. 

2014) . 

It is now well established that there is a strong evolution 
in the normalization of the MS with redshift. The charac¬ 
teristic specific star formation rates (sSFR = SFR/M*) of 
the MS population evolves strongly with redshift, decreas¬ 
ing by a factor of ~ 20 from 2 = 2 to today (e.g., Schreiber 
et al. 2015). All, hydrodynamic simulations of galaxies (e.g., 
Dave et al. 2011; Dekel et al. 2013; Torrey et al. 2014; Sparre 
et al. 2015), semi-analytical models (e.g., Dutton et al. 2010; 
Dave et al. 2012; Mitchell et al. 2014) and analytical models 
(e.g., Bouche et al. 2010; Dekel et al. 2013; Lilly et al. 2013; 
Forbes et al. 2014; Dekel & Mandelker 2014), naturally re¬ 
produce a correlation between SFR and M*. These studies 
show that a natural way to understand the decline of the 
sSFR with time is provided by the predicted decline of gas 
accretion rate onto the galaxies, which itself is closely re¬ 
lated to the evolution of the cosmological specific accretion 
rate into dark matter haloes, which scales as oc (1 -I- z)^'^ at 
a fixed mass in the Einstein-deSitter regime, valid at 2 > 1 
(Neistein et al. 2006; Birnboim et al. 2007; Neistein & Dekel 
2008; Fakhouri & Ma 2009; Genel et al. 2010; Dutton et al. 
2010; Bouche et al. 2010; Tacchella et al. 2013; Lilly et al. 
2013; Dekel et al. 2013). 

As SFGs grow in mass, they seem to propagate along 
the MS, typically not deviating by more than ±0.3 dex from 
the MS ridge^, whose sSFR amplitude steadily declines in 
time. What is the mechanism that keeps the evolving galaxy 
so tightly confined to the vicinity of the MS ridge until it 
quenches and falls below the MS? From the cosmological 
paradigm, dark matter haloes, and hence the central galax¬ 
ies occupying them, form hierarchically - large haloes are 
built from mergers of smaller haloes. One would therefore 
expect that mergers between galaxies would frequently trig¬ 
ger starbursts that would generate larger excursions about 
the MS ridge. However, the small scatter in the MS at multi¬ 
ple redshifts has indicated that most galaxies are not in fact 
experiencing the expected dramatic effects of major merg¬ 
ers (Noeske et al. 2007a,b; Rodighiero et al. 2011), and most 
stars form in “normal” galaxies lying along this relation. The 
SFRs of MS galaxies seem to be sustained for extended peri¬ 
ods of time in a quasi-steady state of gas inflow, gas outflow, 
and gas consumption (Daddi et al. 2010; Bouche et al. 2010; 
Genzel et al. 2010; Tacconi et al. 2010; Dave et al. 2012; Lilly 
et al. 2013; Dekel et al. 2013; Dayal et al. 2013; Feldmann 

2015) . 

Observations indicate that the gas fraction and deple¬ 
tion time tend to vary as a function of sSFR across the MS: 
high sSFR is correlated with high gas fraction and short de- 


^ The ridge is generally the line connecting the medians of the 
sSFR at a given stellar mass, or the points of maximum num¬ 
ber density of galaxies at the given mass. These two definitions 
roughly coincide as the distribution about the ridge is roughly 
symmetric (and log-normal). A detailed definition of the MS ridge 
is provided by Renzini & Peng (2015). 


pletion time (Magdis et al. 2012; Sargent et al. 2014; Huang 
& Kauffmann 2014; Genzel et al. 2015; Silverman et al. 2015; 
Scoville et al. 2015). These gradients may provide a clue for 
understanding the MS width. Therefore, in this paper, we 
focus on galaxy properties as a function of sSFR with respect 
to the MS ridge, rather than the absolute value of sSFR. We 
define the universal MS to be the sSFR with respect to the 
sSFR of the MS ridge. The main questions that we address 
in this paper are: (i) What is the mechanism that confines 
the MS to a small scatter? (ii) What drives the gradients of 
galaxy properties across the universal MS? 

Dutton et al. (2010) used a semi-analytical model for 
disc galaxies to explore the origin of the time evolution and 
scatter of the MS. They find a significant but small scatter 
in their model MS arising from variation in halo concentra¬ 
tion, which in turn causes differences in the mass accretion 
histories between different galaxies of the same halo mass 
(Wechsler et al. 2002). Forbes et al. (2014) presented a toy 
model in which the scatter ultimately arises from the intrin¬ 
sic scatter in the accretion rate, but may be substantially 
reduced depending on the timescale on which the accre¬ 
tion varies compared to the timescale on which the galaxy 
loses gas mass. They show that observational constraints on 
the scatter in galaxy scaling relations can be translated into 
constraints on the galaxy-to-galaxy variation in the outflow 
mass loading factor at fixed mass, and the timescales and 
magnitude of a stochastic component of accretion onto the 
SFGs. 

The key question is which timescale is encoded in 
the MS scatter, i.e., does the MS scatter arise because 
galaxies change their SFR on short timescales (~ 10^ yr), 
intermediate-timescales (~ 10®“® yr), or long timescales 
(~ 10^® yr) (Abramson et al. 2015; Munoz & Peeples 2015). 
If the MS scatter arises due to short term fluctuations in 
the star-formation history, similar mass SFGs mostly grow- 
up together (e.g., Peng et al. 2010; Behroozi et al. 2013). 
On the other hand, if the MS scatter arises due to long term 
fluctuations, similar massive SFGs do not grow-up together 
and key physics lies in what diversifies star-formation histo¬ 
ries (e.g., Gladders et al. 2013; Kelson 2014). 

How SFGs grow their mass during their life on the MS is 
also crucial to understand the build up of the quenched pop¬ 
ulation, and the evolution of its median size. Carollo et al. 
(2013) argue indeed for a straight mass-dependent quench¬ 
ing process (Peng et al. 2010) of M* disc galaxies from the 
MS to the quenched population (with dry mergers playing 
a major role in building the >>M* quenched population, 
whose properties point at a dissipationless process as their 
last step in their assembly histories). A test of this picture is 
to explore how SFGs grow in mass and size on the MS, and 
compare the properties of the massive systems that transi¬ 
tion to the quenched population at the end of their active 
lives. 

We argue here that the confining mechanism and the 
gradients across the MS can be understood in terms of the 
gas regulation, i.e., the balance between inflow rate, SFR 
and outflow rate, of SFGs at high redshift. We emphasise 
the importance of internal physical process, likely driven 
by external events, in addition to global processes (such 
as gas accretion history). Zolotov et al. (2015), analysing 
cosmological zoom-in simulations, have shown that the pro- 
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cesses of gas compaction and subsequent central depletion 
and quenching are frequent in high-z galaxies and are the 
major events in their history. They find that stream-fed, 
highly perturbed, gas-rich discs undergo phases of dissipa¬ 
tive contraction into compact, star-forming systems (“blue 
nuggets”^) at 2 ~ 4 — 2. The compaction is triggered by an 
intense inflow episode, involving mergers, counter-rotating 
streams or recycled gas (Dekel & Burkert 2014), and can be 
associated with violent disc instability (VDI; Noguchi 1999; 
Gammie 2001; Bournaud et al. 2007; Dekel et al. 2009b; 
Burkert et al. 2010; Bournaud et al. 2012; Cacciato et al. 
2012). The peak of gas compaction marks the onset of cen¬ 
tral gas depletion and inside-out quenching. 

Here, we try to learn how this characteristic chain of 
events predicts the gradients across the MS and explains 
the conhnement mechanism. We do this by utilizing the 
same high-resolution, zoom-in, hydro-cosmological, Adap¬ 
tive Mesh Refinement (AMR) simulations as Zolotov et al. 
(2015), of galaxies in the redshift range 2 = 7 to 2 = 1. 
The suite of 26 galaxies analysed here were simulated at a 
maximum resolution of ~ 25 pc including supernova and 
radiative stellar feedback. At 2 ~ 2, the halo masses are in 
the range Mvu ~ 10^^“^^ Mq and the stellar masses are 
in the range M* ~ 20®'^“^°'® Mq. With these simulations, 
we focus on the global physical properties of galaxies on the 
MS. 

This paper is organized as follows. In Section 2, we give 
a brief overview of the simulations. In Section 3, we investi¬ 
gate the gas content of the simulated galaxies. In Section 4, 
we define the MS, and in Section 5, we determine galaxy 
properties across the MS. The core of this paper is Section 6, 
where we explain the confinement mechanism of the MS. We 
discuss implications from our MS paper for the cessation of 
star formation in galaxies and we highlight several caveats 
of our analysis in Section 7. We summarize our results in 
Section 8. 


2 SIMULATIONS 

We use zoom-in hydro-cosmological simulations of 26 mod¬ 
erately massive galaxies, a subset of the 35-galaxy VELA sim¬ 
ulation suite. The details the VELA simulations are presented 
in Ceverino et al. (2014) and Zolotov et al. (2015). Zolotov 
et al. (2015) and Tacchella et al. (2015a) used the same 
sample of 26 simulations and investigated similar questions 
concerning compaction and quenching. Zolotov et al. (2015) 
focused on the evolution of the global properties of the galax¬ 
ies and their cores as they go through the compaction and 
quenching phases, and Tacchella et al. (2015a) addresses the 
evolution of surface density profile of these galaxies during 
these phases. Additional analysis of the same suite of sim¬ 
ulations are discussed in Moody et al. (2014), Snyder et al. 
(2015) and Ceverino et al. (2015b). In this section, we give 
an overview of the key aspects of the simulations. 


^ In this paper, we refer to compact, SFGs as a blue nuggets. 
Such galaxies have a high density in their cores, both in stellar 
mass and gas density. Note that blue nuggets could actually be 
quiet red due to dust. 


2.1 Cosmological Simulations 

The VELA simulations utilize the Adaptive Refinement Tree 
(ART) code (Kravtsov et al. 1997; Kravtsov 2003; Ceverino 
& Klypin 2009), which accurately follows the evolution of 
a gravitating A-body system and the Eulerian gas dynam¬ 
ics. All the simulations were evolved to redshifts 2 < 2, 
and several of them were evolved to redshift 2 = 1, with 
an AMR maximum resolution of 17 — 35 pc at all times, 
which is achieved at densities of ~ 10“^ — 10® cm“®. In 
the circumgalactic medium (at the virial radius of the dark- 
matter halo), the median resolution amounts to ~ 500 pc. 
Beyond gravity and hydrodynamics, the code incorporates 
the physics of gas and metal cooling, UV-background pho¬ 
toionization, stochastic star formation, gas recycling and 
metal enrichment, and thermal feedback from supernovae 
(Ceverino et al. 2010, 2012), plus a new implementation of 
feedback from radiation pressure (Ceverino et al. 2014). 

We use the CLOUDY code (Ferland et al. 1998) to cal¬ 
culate the cooling and heating rates for a given gas density, 
temperature, metallicity, and UV background, assuming a 
slab of thickness 1 kpc. We assume a uniform UV back¬ 
ground, following the redshift-dependent Haardt & Madau 
(1996) model. An exception is at gas densities higher than 
0.1 cm“®. At these densities, we use a substantially sup¬ 
pressed UV background (5.9 x 10® erg s“® cm“® Hz“®) in 
order to mimic the partial self-shielding of dense gas, allow¬ 
ing dense gas to cool down to temperatures of ~ 300 K. The 
equation of state is assumed to be that of an ideal mono- 
atomic gas. Artificial fragmentation on the cell size is pre¬ 
vented by introducing a pressure floor, which ensures that 
the Jeans scale is resolved by at least 7 cells (see Ceverino 
et al. 2010). 

We assume that star formation occurs at densities above 
a threshold of 1 cm“® and at temperatures below lO'^ K. 
Most stars (> 90 %) form at temperatures well below 10® K, 
and more than half of the stars form at 300 K in cells where 
the gas density is higher than 10 cm“®. We use a stochas¬ 
tic star-formation model, where star formation occurs in 
timesteps of dtsp = 5 Myr. The probability to form a stellar 
particle in a given timestep is 

P = min ( 0 . 2 ,, (1) 

V V 1000 cm-3 J ^ ’ 

The single stellar particle has a mass equal to 

m* = rrigas ~ 0.42mgas (2) 

T 

where rrigas is the mass of gas in the cell where the particle is 
being formed and t is 12 Myr. We assume a Chabrier (2003) 
initial mass function. This stochastic star-formation model 
yields a star-formation efficiency per free-fall time of ~ 2 %. 
At the given resolution, this efficiency roughly mimics the 
empirical Kennicutt-Schmidt law (Kennicutt 1998). As a re¬ 
sult of the universal local SFR law adopted, the global SFR 
follows the global gas mass (see Figs. 2 and 3 in Zolotov et al. 
2015). Observationally, a universal, local SFR law in which 
the star formation rate is simply ~ 1 % of the molecular 
gas mass per local free-fall time fits galactic clouds, nearby 
galaxies, and high-redshift galaxies (Krumholz et al. 2012). 
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The suite of 26 simulated galaxies. 


Galaxy 

Mvir 

10^2 Mq 
{ z = 2) 

M* 

10l° Mq 
( z = 2) 

Afgas 
10l° Mq 
iz = 2) 

SFR 
Mq/jv 
{ z = 2) 

sSFR 
Gyr-i 
(z = 2) 

-^vir 
kpc 
(^ = 2) 

Rm 

kpc 
{z = 2) 

^fin 

•2^ fin 

01 

0.16 

0.22 

0.12 

2.65 

1.20 

58.25 

1.06 

0.50 

1.00 

02 

0.13 

0.19 

0.16 

1.84 

0.94 

54.50 

2.19 

0.50 

1.00 

03 

0.14 

0.43 

0.10 

3.76 

0.87 

55.50 

1.7 

0.50 

1.00 

06 

0.55 

2.22 

0.33 

20.72 

0.93 

88.25 

1.06 

0.37 

1.70 

07 

0.90 

6.37 

1.42 

26.75 

0.42 

104.25 

2.78 

0.50 

1.00 

08 

0.28 

0.36 

0.19 

5.76 

1.58 

70.50 

0.76 

0.50 

1.00 

09 

0.27 

1.07 

0.31 

3.97 

0.37 

70.50 

1.82 

0.39 

1.56 

10 

0.13 

0.64 

0.11 

3.27 

0.51 

55.25 

0.53 

0.50 

1.00 

11 

0.27 

1.02 

0.58 

17.33 

1.69 

69.50 

2.98 

0.46 

1.17 

12 

0.27 

2.06 

0.19 

2.91 

0.14 

69.50 

1.22 

0.39 

1.56 

13 

0.31 

0.96 

0.98 

21.23 

2.21 

72.50 

3.21 

0.39 

1.56 

14 

0.36 

1.40 

0.59 

27.61 

1.97 

76.50 

0.35 

0.42 

1.38 

15 

0.12 

0.56 

0.14 

1.71 

0.30 

53.25 

1.31 

0.50 

1.00 

20 

0.53 

3.92 

0.48 

7.27 

0.19 

87.50 

1.81 

0.44 

1.27 

21 

0.62 

4.28 

0.57 

9.76 

0.23 

92.25 

1.76 

0.50 

1.00 

22 

0.49 

4.57 

0.21 

12.05 

0.26 

85.50 

1.32 

0.50 

1.00 

23 

0.15 

0.84 

0.19 

3.32 

0.39 

57.00 

1.38 

0.50 

1.00 

24 

0.28 

0.95 

0.28 

4.39 

0.46 

70.25 

1.79 

0.48 

1.08 

25 

0.22 

0.76 

0.08 

2.35 

0.31 

65.00 

0.82 

0.50 

1.00 

26 

0.36 

1.63 

0.25 

9.76 

0.60 

76.75 

0.76 

0.50 

1.00 

27 

0.33 

0.90 

0.52 

8.75 

0.97 

75.50 

2.45 

0.50 

1.00 

29 

0.52 

2.67 

0.39 

18.74 

0.70 

89.25 

1.96 

0.50 

1.00 

30 

0.31 

1.71 

0.41 

3.84 

0.22 

73.25 

1.56 

0.34 

1.94 

32 

0.59 

2.74 

0.37 

15.04 

0.55 

90.50 

2.6 

0.33 

2.03 

33 

0.83 

5.17 

0.45 

33.01 

0.64 

101.25 

1.22 

0.39 

1.56 

34 

0.52 

1.73 

0.42 

14.79 

0.85 

86.50 

1.9 

0.35 

1.86 


Table 1. Quoted are the total virial mass, Mvir, the stellar mass, M*, the gas mass, Mgas, the star formation rate, SFR, the specific star 
formation rate, sSFR, the virial radius, Rvin the effective stellar (half-mass) radius, -Rm> a-U a,t z = 2, and the final simulation snapshot, 
Ufin, and redshift, za^. The M*, Mg, SFR, and sSFR are measured within a radius of 0.2 X Rvir- 


The thermal stellar feedback model releases energy from 
stellar winds and snpernova explosions as a constant heating 
rate over 40 Myr following star formation. The heating rate 
due to feedback may or may not overcome the cooling rate, 
depending on the gas conditions in the star-forming regions 
(Dekel & Silk 1986; Ceverino & Klypin 2009). Note that no 
artificial shutdown of cooling is implemented in these simu¬ 
lations. The effect of runaway stars is included by applying 
a velocity kick of ~ 10 km s“^ to 30 % of the newly formed 
stellar particles. The code also includes the later effects of 
Type la supernova and stellar mass loss, and it follows the 
metal enrichment of the ISM. 

Radiation pressure is incorporated through the addition 
of a non-thermal pressure term to the total gas pressure in 
regions where ionizing photons from massive stars are pro¬ 
duced and may be trapped. This ionizing radiation injects 
momentum around massive stars, pressurizing star-forming 
regions, as described in Appendix B of Agertz et al. (2013). 
We assume an isotropic radiation field within a given cell and 
that the radiation pressure is proportional T • m*, where m^, 
is the mass of stars and T is the luminosity of ionizing pho¬ 
tons per unit stellar mass. The value of T is taken from the 
stellar population synthesis code, STARBURST99 (Leitherer 
et al. 1999). We use a value of T = 10^® erg s“^ which 

corresponds to the time-averaged luminosity per unit mass 
of the ionizing radiation during the first 5 Myr of evolution 
of a single stellar population. After 5 Myr, the number high 


mass stars and ionizing photons declines significantly. Fur¬ 
thermore, the significance of radiation pressure also depends 
on the optical depth of the gas within a cell. We use a hy¬ 
drogen column density threshold, N = 10^^ cm“^, above 
which ionizing radiation is effectively trapped and radiation 
pressure is added to the total gas pressure. This value corre¬ 
sponds to the typical column density of cold neutral clouds, 
which host optically-thick column densities of neutral hy¬ 
drogen (Thompson et al. 2005). Summarizing, our current 
implementation of radiation pressure adds radiation pres¬ 
sure to the total gas pressure in the cells (and their clos¬ 
est neighbours) that contain stellar particles younger than 
5 Myr and whose column density exceeds cm“^. 

2.2 Limitation of the current Simulations 

The cosmological simulations used in this paper are state- 
of-the-art in terms of high-resolution AMR hydrodynamics 
and the treatment of key physical processes at the subgrid 
level, highlighted above. Specifically, these simulations trace 
the cosmological streams that feed galaxies at high redshift, 
including mergers and smooth flows, and they resolve the 
VDI that governs high-z disc evolution and bulge formation 
(Ceverino et al. 2010, 2012, 2015a; Mandelker et al. 2014). 

Like other simulations, the current simulations are not 
yet doing the best possible job treating the star formation 
and feedback processes. As mentioned above, the code as- 
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sumes a SFR efficiency per free fall time that is rather re¬ 
alistic, but it does not yet follow in detail the formation of 
molecules and the effect of metallicity on SFR (Krumholz 
& Dekel 2012). Additionally, the resolution does not allow 
the capture of Sedov-Taylor adiabatic phase of supernova 
feedback. The radiative stellar feedback assumed no infrared 
trapping, in the spirit of low trapping advocated by Dekel & 
Krumholz (2013) based on Krumholz & Thompson (2012). 
Other works assume more significant trapping (Murray et al. 
2010; Krumholz & Dekel 2010; Hopkins et al. 2012), which 
makes the assumed strength of the radiative stellar feed¬ 
back here lower than in other simulations. Finally, AGN 
feedback and feedback associated with cosmic rays and mag¬ 
netic fields are not yet incorporated. Nevertheless, as shown 
in Ceverino et al. (2014), the star formation rates, gas frac¬ 
tions, and stellar-to-halo mass fractions are all in the ball¬ 
park of the estimates deduced from abundance matching, 
providing a better match to observations than earlier simu¬ 
lations. 

The uncertainties and any possible remaining mis¬ 
matches between simulation and observations by a factor 
of order 2 are comparable to the observational uncertain¬ 
ties. For example, the stellar-to-halo mass fraction is not 
well constrained observationally at 2 ~ 2. Recent estimates 
by Burkert et al. (2015) (see their Fig. 5) based on the ob¬ 
served kinematics of 2 ~ 0.6 — 2.8 SFGs reveal significantly 
larger ratios than the estimates based on abundance match¬ 
ing (Conroy & Wechsler 2009; Moster et al. 2010, 2013; 
Behroozi et al. 2010, 2013) at Mvu < 10^^ Mq. In Sec¬ 
tion A of the appendix, we present a detailed comparison 
of the stellar-to-halo mass relation of our simulations and 
the observational data (Figure Al). We conclude that our 
simulations produce stellar-to-halo mass ratios that are in 
the ballpark of the values estimated from observations, and 
within the observational uncertainties, thus, indicating that 
our feedback model is adequate. 

It seems that in the current simulations, the compaction 
and the subsequent onset of quenching occur at cosmological 
times that are consistent with observations (see Fig. 12 of 
Zolotov et al. 2015 and Fig. 2 of Barro et al. 2013). With 
some of the feedback mechanisms not yet incorporated (e.g., 
AGN feedback), full quenching to very low sSFR values may 
not be fully reached in many galaxies by the end of the 
simulations at 2 ~ 1. In this work, we adopt the hypothesis 
that the simulations grasp the qualitative features of the 
main physical processes that govern galaxy evolution. 


2.3 Galaxy Sample and Properties 

The initial conditions for the simulations are based on dark- 
matter haloes that were drawn from dissipationless N-body 
simulations at lower resolution in three comoving cosmolog¬ 
ical boxes (box-sizes of 10, 20, and 40 Mpc/h). The assumed 
cosmology is the standard ACDM model with the WMAP5 
values of the cosmological parameters, namely = 0.27, 
Ha = 0.73, Hi, = 0.045, h = 0.7 and cts = 0.82 (Komatsu 
et al. 2009). Each halo was selected to have a given virial 
mass at 2 = 1 and no ongoing major merger at 2 = 1. 
This latter criterion eliminates less than 10 % of the haloes, 
which tend to be in a dense proto-cluster environment at 
2 ~ 1. The target virial masses at 2 = 1 were selected in the 


range Mvir = 2 x 10^^ — 2 x 10^^ Mq, about a median of 
4.6 X 10^^ Mq. If left in isolation, the median mass at 2 = 0 
was intended to be ~ 10^^ Mq. In practice, the actual mass 
range is broader, with some of the haloes merging into more 
massive haloes that eventually host groups at 2 = 0. 

From the suite of 35 galaxies, we have excluded three 
low-mass galaxies that have a quenching attempt that brings 
them considerably below the MS (> 0.8 dex) for a short 
period before they recover back to the MS. This is probably 
a feature limited to very low mass galaxies, and we do not 
address it any further here. Furthermore, we have excluded 
six galaxies which have not been simulated down to 2 = 2.0. 
Therefore, our final sample consists of 26 galaxies. The virial 
and stellar properties are listed in Table 1. This includes 
the total virial mass Mvir, the galaxy stellar mass M*, the 
gas mass Mgas, the SFR, the sSFR, the virial radius i?vir, 
and the effective, and half-mass radius Rm, all at 2 = 2. 
The latest time of analysis for each galaxy in terms of the 
expansion factor, agn, or redshift, 2 fin, is provided. 

The virial mass Mvir is the total mass within a sphere 
of radius Rvir that encompasses an overdensity of A( 2 ) = 
(187r^ — 82 Ha( 2:) — 39HA(2)^)/Hm(2:), where Ha( 2) and 
H„i( 2 ) are the cosmological parameters at 2 (Bryan & Nor¬ 
man 1998; Dekel & Birnboim 2006). The stellar mass of the 
galaxy, M*, is the instantaneous mass in stars (after the 
appropriate mass loss), measured within a sphere of radius 
0.2 X Rvir about the galaxy centre. The gas mass, Mgas, is 
the cold gas mass within the same sphere, i.e., the gas with 
a temperature below 10"* K. Measuring these global quanti¬ 
ties within different radii (0.1 — 0.3 x Rvir or even fixing it 
to 10 kpc) does not change the main findings of this paper. 
Throughout this paper, all the quoted gas properties refer 
to the cold gas component. The effective radius Rm is the 
three-dimensional half-mass radius corresponding to M*. 

The SFR is obtained by SFR = (M*(fage < 
tmax)/tmax)t„ax> where M*(tage < tmax) is the mass in stars 
younger than tmax within a sphere of radius 0.2 x Rvir about 
the galaxy centre. The average (•)t„aj, is obtained for tmax 
in the interval [40, 80] Myr in steps of 0.2 Myr in order to 
reduce fluctuations due to a ~ 5 Myr discreteness in stellar 
birth times in the simulation. The tmax in this range are long 
enough to ensure good statistics. 

We start the analysis at the cosmological time corre¬ 
sponding to expansion factor a = 0.125 (redshift 2 = 7). 
At earlier times, the fixed resolution scale typically corre¬ 
sponds to a larger fraction of the galaxy size, so the detailed 
galaxy properties may be less accurate. As visible in Ta¬ 
ble 1, most galaxies reach a = 0.50 (2 = 1). Each galaxy 
is analysed at output times separated by a constant inter¬ 
val in a, Aa = 0.01, corresponding at 2 = 2 to ~ 100 Myr 
(roughly half the orbital time at the disc edge). For six galax¬ 
ies, namely 11, 12, 14, 25, 26, and 27, ~ 20-times higher res¬ 
olution snapshots are available, for which the output times 
separated by Ao = 0.0005 — 0.0007. In Appendix B, we show 
that our standard snapshots are tracing the main evolution¬ 
ary pattern, and the high temporal resolution snapshots con¬ 
firm our main findings. 
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Figure 1. Gas fraction versus sSFR. The circles correspond to 
all the snapshots of all the galaxies in the simulated sample. The 
colour coding corresponds to redshift according to the colour bar. 
We see that /gas and sSFR are systematically decreasing with 
cosmic time. The simulations are compared to the observations 
of Tacconi et al. (2013), marked by stars. The errorbar indicates 
the uncertainty in the observations. The simulated and observed 
galaxies span a similar locus in this plane, despite the fact that the 
simulations somewhat underestimate both the gas fraction and 
the sSFR (see text). The black solid line indicates the best-fitting 
depletion time at z 2.5 of the simulations (tuep = 480ibl0 Myr), 
which is shorter than the observational estimate of tjep = 700 it 
200 Myr shown as a black dashed line. 

3 GAS CONTENT IN SIMULATED GALAXIES 

As mentioned before, the star formation in the simulations 
is driven by the content of dense and cold gas in the galaxies. 
It is therefore important to investigate the gas properties of 
our simulated galaxies, before analysing the evolution and 
shape of the MS. Figure 1 shows the gas fraction, /gas = 
Mgas/(Mgas + A/*), as a function of sSFR. The individual 
points refer to the individual snapshots, ranging from 2 = 6 
to 1 as indicated by the colour coding. The high redshift 
(2 > 3) galaxies of typical masses below ~ 10®'® Mq and 
sSFR ~ 2 Gyr“^ tend to have a high gas fraction, /gas > 0.4. 
Towards lower redshifts (2 ~ 1 — 3), when the galaxies in our 
sample also grow to higher masses, the /gas and the sSFR 
continuously decline with increasing cosmic time. At 2 ~ 2, 
our simulated 10^® Mq galaxies have /gas = 0.20to',o9 (see 
below). 

The quantities /gas and sSFR are related through 


Mgas + M* 1 -b (tdep X sSFR)-i ’ ^ ^ 

where tdep = Mgas/SFR is the depletion time®. From Fig- 

® The depletion time (sometimes also called the gas consumption 
time) can be written as tdep = Mgas/SFR = tff/Eff, where Sfj is 
the SFR efficiency and tfj is the free-fall time in the star-forming 
regions. The SFR efficiency can be argued to be constant in all 
star-forming environments, eg ~ 0.01, such that the variation in 
tdep mostly reflects variations in fg (Krumholz et al. 2012). 


ure 1 , we can estimate a global tdep by fitting all our sim¬ 
ulated galaxies at all redshifts together. We find tdep = 
400 ± 10 Myr in the simulations. Considering only galax¬ 
ies at 2 2.5, we find tdep = 480 ± 10 Myr. As a by-product 

to our current analysis, from the fact that this timescale is 
much shorter than the Hubble time for all SFGs between 
2 = 1 — 4, we can conclude that fresh gas must be supplied 
with a fairly high duty cycle over several billion years. 

We compare the simulations with measurements of Tac¬ 
coni et al. (2013) who present a CO 3 — 2 survey of molec¬ 
ular gas properties in massive galaxies at 2 ~ 1 — 3. They 
provided 52 CO detections in galaxies with logj^Q M*/Mq > 
10.4 and logj^p SFR/(M 0 yr“®) > 1.5. The galaxies were se¬ 
lected to sample the complete range of star formation rates 
within the aforementioned stellar mass limit, ensuring a rel¬ 
atively uniform sampling of the MS. They infer average gas 
fractions of /gas ~ 0.33 at 2 ~ 1.2 and ~ 0.47 at 2 ~ 2.2 for 
the given masses, and an overall drop in /gas with M* at a 
given redshift. 

Both the observed and the simulated galaxy sample 
are incomplete and the comparison between the two has to 
be interpreted with caution. Comparing the Tacconi et al. 
(2013) measurements with our simulations, we find that 
galaxies in the simulations have an average gas deficit at 
a given stellar mass of 20% — 40% in comparison with 
the observations. However, at a given sSFR, we find good 
agreement between observations and simulations (Figure 1). 
Tacconi et al. (2013) inferred an average depletion time 
tdep = 700 ± 200 Myr at 2 = 1 — 3, which may be slightly 
larger than our estimate of ~ 480 Myr from the simulations. 
The good agreement is found because at a given M*, the 
SFR in the simulations is also underestimated by a similar 
multiplicative factor as the gas density. We conclude that al¬ 
though the simulated galaxies may evolve a bit ahead of cos¬ 
mic time, they nevertheless reproduce the qualitative trends 
seen in the observations. 


4 THE STAR-FORMING MAIN-SEQUENGE 

In this section, we identify the MS in the simulations and 
determine the strong systematic time evolution of its ridge 
and its weak mass dependence, sSFRms(A/*, 2 ). We are mo¬ 
tivated by the hypothesis that the sSFRMs(Af*, 2 ) roughly 
follows the average specific accretion rate of mass into dark- 
matter haloes, and test this hypothesis. After obtaining the 
best-fitting MS in the simulations, we compare it with ob¬ 
servations. 

4.1 Definition of the MS in the Simulations 

The average specific accretion rate of mass into haloes of 
mass Mh at 2 can be approximated by an expression of the 
form: 

^:^sg.Mf 2 .(l + 2 )^ (4) 

where M 12 = Mg/IO^^M©. In the ACDM cosmology in the 
Einstein-deSitter regime (valid at 2 > 1), simple theoreti¬ 
cal arguments show that p —>■ 5/2 (see Dekel et al. 2013). 
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Figure 2. Evolution of the MS in the simulations and the universal MS. The colour refers to redshift according to the colour bar. 
Left panel: sSFR as a function of stellar mass M*. The individual points show the galaxies from the simulations at different redshifts 
{z = G — 1). The solid lines show the evolution of the MS ridge sSFRms (Equation 5). Right panel: The universal MS, i.e., the distance 
from the MS ridge, Ams = logjQ(sSFR/sSFR]y[s), as a function of M* and z. The solid line marks the MS ridge and the gray shaded 
region indicates the ±0.3 dex scatter of the MS. The simulations recover the confinement to a narrow MS, with a bending down at the 
massive end at late times. 


Neistein & Dekel (2008) showed that, for the ACDM power- 
spectrum slope on galactic scales, P is small, /3 ~ 0.14. An 
appropriate value for the normalization factor Sh provides a 
match better than 5% at 2 > 1 to the average evolution in 
cosmological N-body simulations (Dekel et al. 2013). 

To constrain the evolution of the MS ridge with cosmic 
time, we adopt for the sSFR the same functional form as in 
Equation 4, 

sSFRms(M*, z) = Sb • j • (1 + C" Gyr"\ (5) 

The three free parameters (sb, P, and fi) are not independent 
of each other. In particular, the time evolution of galaxies 
is characterized by both /3 (stellar mass dependence) and 
pL (redshift evolution). We therefore assume P = 0.14 and 
pi = 5/2, i.e., the values of the specific halo mass accretion 
rate. This choice of P and pi is motivated by our data. Fitting 
P in bins of redshift (bins of 0.5 in the range 2 = 6 — 3), we 
find as best-fitting P = 0.12 ± 0.06. Using P = 0.12 and 
htting pi, we find pi = 2.5 ± 0.6. 

Since we adopt P = 0.14 and pi — 5/2, the only free pa¬ 
rameter in the MS-scaling is the overall normalisation which 
we parametrise with Sb. We determine Sb by a least-square 
ht to logj^Q sSFR using all the snapshots in the redshift range 
6 — 3 for all the galaxies in our sample, excluding all galax¬ 
ies that are below the 16%-percentile or above the 84%- 
percentile at a given redshift, thus focusing on the inner 
±1(T about the median. The choice of this redshift range is 
motivated by the fact that all of our galaxies in the range 
2 = 6 — 3 are star-forming. We find Sb = 0.046 ± 0.002. Al¬ 
ternatively, if we replace the redshift cut by a cut in stellar 


surface density within the inner 1 kpc, SM*,ikpc < 10® M©, 
as in Zolotov et al. (2015), we obtain a value of Sb that 
is consistent within the uncertainty with the above value, 
which we adopt here. 

Figure 2 shows the evolution of the best-htting MS ridge 
in the simulations. The left panel shows the MS ridge with 
solid lines at redshifts 2 = [1.0, 2.0,3.0,4.0, 5.0, 6.0], indi¬ 
cated by the different colours following the colour bar. The 
individual points refer to all the snapshots of all the sample 
galaxies at redshifts 1 ^ 2 ^ 6. The right panel of Figure 2 
shows the universal MS, namely all snapshots after scaling 
the sSFR according to the scaling of the MS ridge. The hg- 
ure thus shows the distance of each galaxy from the MS 
ridge, Ams = logig (sSFR/sSFRms). The grey shaded area 
indicates a scatter of ±0.3 dex, indicative of the scatter of 
the MS. 

Two low-mass galaxies become sub-MS at 2 ~ 4 —3 and 
then return to the MS (galaxies 10 and 24). These galaxies 
live through a short (< 200 Myr) phase of nearly no accre¬ 
tion of gas that reduces their SFR significantly. This quench¬ 
ing attempt is terminated by a sudden accretion of fresh 
gas, usually initiated by a merger. This brings the galaxy 
back into the MS within less than 400 Myr. Such episodes 
of quenching attempts seem to occur more frequently in 
low-mass galaxies at high- 2 , but they are rare above the 
mass threshold adopted for the sample analysed in this pa¬ 
per (Section 2.3). 

At later epochs (2 < 2), for several galaxies, the quench¬ 
ing process is continuous over several 100 Myr up to several 
Gyr, indicating a successful quenching as opposed to a short¬ 
term quenching attempt at high- 2 . Only 3 galaxies fall 1 dex 
below the MS by 2 = 1, i.e. quench their star-formation 


MNRAS 000, 1-27 (2016) 








S. Tacchella, A. Dekel, C. M. Carollo, et al. 



redshift redshift 


Figure 3. Comparison of the sSFR in the MS ridge from the 
simulations with observations for galaxies at two mass bins 
(logj^QM*/M 0 ~ 10.5 in the left and extrapolated to 11.5 in 
the right panel). The red solid line shows the MS ridge from our 
simulations, i.e., the best-fitting sSFRms (Equation 5). The blue 
dashed-dotted, green dashed, and purple dotted lines show the 
best-fitting MS ridges of the observations by Lilly et al. (2013), 
Speagle et al. (2014), and Schreiber et al. (2015), respectively. 
The shade regions enclosing the lines of the observations with a 
0.3 dex width show the approximate observational uncertainties. 
At 2 ~ 2, the simulations match the observations in the high- 
mass end, but they are an underestimate by 0.3 dex for the lower 
masses. 


substantially. When these galaxies fall below the MS, their 
stellar mass roughly stays constant. As mentioned in Sec¬ 
tion 2.2, full quenching to very low sSFR is not reached 
because the feedback may still be underestimated (e.g., the 
adiabatic phase of supernova feedback is not resolved, and 
AGN feedback is not yet incorporated). However, 12 out of 
our 26 galaxies drop to more than 1 cr below the MS ridge by 
2 = 1, and all are continuously moving downward in Ams for 
the last several Gyr, indicating that they are in a long-term 
quenching process. 


4.2 Comparison of the MS in the Simulations 
with Observations 

Figure 3 compares the sSFR amplitude of the MS ridge 
from the simulations with the one deduced from observations 
(Lilly et al. 2013; Speagle et al. 2014; Schreiber et al. 2015) 
for two different stellar masses. Our simulations, the red 
curves, show the best-fitting MS introduced above (Equa¬ 
tion 5), extrapolated to masses and redshifts beyond what 
we actually have in the simulations (e.g., we do not have 
any 10^^ Mq galaxies at 2 > 2 in the simulations). The Lilly 
et al. (2013) line is based on best-fitting to the cumulated 
data of Noeske et al. (2007a); Elbaz et al. (2007); Daddi 
et al. (2007); Pannella et al. (2009) and Stark et al. (2013). 
The line of Speagle et al. (2014) uses a compilation of 25 
observational studies from the literature out to 2 ~ 6. After 
converting all observations to a common set of calibrations, 
they find a remarkable consensus among MS ridge observa¬ 
tions, with ~ 0.1 dex for the 1 a inter-publication scatter. 
Schreiber et al. (2015) presented an analysis of the deep¬ 
est Herschel images obtained within the GOODS-Herschel 
and CANDELS-Herschel key programs. This allowed them 
to measure SFR based on direct ultraviolet and far-infrared 


ultraviolet-reprocessed light and to determine the evolution 
of the MS at 2 = 0 — 4. 

In Figure 3, at high masses (right panel), the MS from 
observations and our simulations are consistent, especially 
at intermediate redshifts, 2 ~ 1 —3. The differences at higher 
redshifts may reach the levels of ~ 0.3 dex, i.e., of the same 
level as potential systemic errors in both SFR and M* de¬ 
duced from observations, and comparable to the difference 
between the different compilations of observations. At lower 
masses (left panel), the MS of the simulations is in less good 
agreement with observations, especially at 2 = 1 — 3, where 
the simulations lie 0.2 — 0.4 dex below the observational es¬ 
timates. 

The larger difference towards lower stellar masses be¬ 
tween the MS of the simulations and the observations can 
be partly explained by a difference in the logarithmic slope 
of the mass dependence: /3 = 0.14 for our simulations (and 
the dark matter haloes, e.g., Neistein & Dekel 2008), where 
/3obs ~ —0.1 to 0.0 in the observations (e.g., Schreiber et al. 
2015). This small difference may be relevant for some physi¬ 
cal processes (e.g., Bouche et al. 2010), but for our purpose, 
given the rather limited mass range in our simulated sample, 
it does not make a difference for the nature of the evolution 
about the MS. Despite these differences between observed 
and simulated MS, we hypothesize that a qualitative study 
of the physical processes governing the evolution of galaxies 
with respect to the MS ridge can be meaningful once treated 
self-consistently using the MS as defined in the simulations 
above. 

As noted before, the cosmological rates of evolution 
of the average observed and simulated sSFR are consistent 
with the average specific accretion rate of mass into haloes 
as expressed in Equation 4 (Bouche et al. 2010; Dekel et al. 
2013; Lilly et al. 2013). However, there are certain differences 
between these two quantities. The sSFR (with Sb ~ 0.046 
Gyr“^) appears to be somewhat higher than the specific 
halo mass accretion rate (sAR, with Sh ~ 0.03 Gyr“^, as 
measured from another set of simulations in Dekel et al. 
2013) over a wide range of redshifts, indicating a factor of 
~ 1.5 shorter e-folding time for the build up of stars com¬ 
pared with that of the dark matter haloes. This is consis¬ 
tent with predictions from a bathtub toy model by Lilly 
et al. (2013) and Dekel & Mandelker (2014) that predict 
sSFR ~ l//star • sAR ~ 1.5 • sAR at high 2 , where /star is 
the mass fraction retained in long-lived stars after stellar 
mass loss. Furthermore, Lilly et al. (2013) pointed out that 
the difference in /3 between the sAR and the sSFR could 
be a result of the curvature in the /star(Af*) relation that 
can itself be traced to the curvature in the mass-metallicity 
relation. Nevertheless, the mass dependence is rather weak 
either way, and it has only a secondary effect on our current 
study. 

4.3 Scatter of the MS 

Observations typically reveal a MS scatter of ums — 0.3 dex 
(e.g., Noeske et al. 2007a; Rodighiero et al. 2011; Whitaker 
et al. 2012; Schreiber et al. 2015). Glearly, the measured 
scatter depends on the exact selection criteria for the SFGs 
as well as on the uncertainties in the stellar mass and SFR 
indicators. 
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As discussed by Speagle et al. (2014), each MS ob¬ 
servation is measured within a predefined redshift window. 
Because of this, the observed scatter about the MS, (tms, 
is actually the intrinsic scatter about the MS convolved 
with the MS’s evolution within the given time interval. The 
(Tms measured from the observations are therefore overesti¬ 
mates. Speagle et al. (2014) deconvolved the observed scat¬ 
ter by using their best fits to approximate the MS evolu¬ 
tion within each time interval, and subtracted this evolution 
from the observed scatter. Furthermore, the observation- 
induced scatter is also taken into account. They find that 
the true, observation-corrected scatter is estimated to be 
~ 0.2 — 0.3 dex, respectively. They find the scatter to be 
roughly constant with cosmic time. 

In the simulations, we can directly measure the true 
scatter about the MS. We measure a scatter of (Jms = 
0.27 ± 0.01 dex for 2 = 3 — 6 (error obtained from boot¬ 
strapping), i.e., consistent with the observational estimates. 
We find a weak trend with cosmic time: the scatter increases 
from 0.25^0 01 dex at 2 ~ 5 to 0.31 ±0.01 dex at 2 ~ 3. This 
may reflect more contamination by quenching galaxies at 
lower redshifts, leading to a bend of the MS downwards at 
the massive end. 


5 GALAXY PROPERTIES ACROSS THE MS 

In this section, we measure galaxy properties in the simula¬ 
tions as they evolve along and across the MS. We attempt 
to correlate the galaxy migration above and below the MS 
ridge with the major events of compaction, depletion, and 
quenching during the galaxy’s evolution. We first look at 
a few galaxies individually. Afterwards, we determine gra¬ 
dients across the MS from all galaxies in the simulations 
and compare them with recent observations by Genzel et al. 
(2015). Finally, we investigate the timescale for the oscilla¬ 
tion along the MS ridge by conducting a Fourier analysis. 

5.1 Individual Galaxies 

In Figures 4, 5, and 6 , we show the evolution of gas density 
within the central 1 kpc (pgas.i), the overall gas to stellar 
mass ratio within the galaxy (/gs = Mgas/A4*), and the de¬ 
pletion time (tdep; Equation 3), for eight galaxies along the 
MS. In the three figures, we plot the distance from the MS, 
Ams = logio (sSFR/sSFRms), as a function of stellar mass, 
M*. The grey vertical line indicates the stellar mass when 
the galaxy’s halo reaches Mvir = ® Mq, argued to mark 

the threshold halo mass for virial shock heating and thus 
long-term quenching (Dekel & Birnboim 2006; Zolotov et al. 
2015). 

Since we wish to quantify the gradients of Pgas.i, /gs, 
and tdep across the universal MS, we have to correct for their 
intrinsic cosmic time evolution. For example, we saw in Sec¬ 
tion 3 that the gas content of a galaxy is a strong function 
of redshift and stellar mass. Since we follow a galaxy sample 
through cosmic time, the time evolution is naturally associ¬ 
ated with growth of stellar mass. We quantify this evolution 
with cosmic time and stellar mass in detail in Section 5.2 
and in Appendix C. Briefly, we determine the average evo¬ 
lution by the functions f{z) = ^ ± x log]^Q(l ± 2 ) and 


g{MP) = 77 ± 7 X (logj^Q(M*) — 10.5) deduced for the quan¬ 
tity of interest from all galaxies that lie close to the MS 
ridge (|Ams| < 0.15), assuming that f{z) and g{M,,) are in¬ 
dependent from each other. The best-fitting parameters can 
be found in Table 2. We then divide all individual measure¬ 
ments by f{z) and g{M,,). 

The eight example galaxies shown in Figures 4, 5, and 
6 are prototypical for two mass bins at 2 = 2 : the upper 
four galaxies (11, 14, 25, 27) belong to the low-mass bin 
(logjg M^/Mq < 10.2), while the lower four (07, 12, 26, 
29) are from the high-mass bin (logj^g M*/M 0 > 10.2). We 
caution that the galaxies which are in a certain mass bin at 
2 = 2 are not necessarily in the same mass bin at all times. 
We make this division in order to differentiate more evolved 
versus less evolved galaxies at a certain epoch, motivated 
by the finding of Zolotov et al. (2015) that more massive 
galaxies at 2 = 2 tend to quench earlier and more decisively. 

The first insight from these plots is that the individual 
galaxies tend not to be super-MS or sub-MS at all times; 
they evolve through super-MS and sub-MS phases. Most 
galaxies actually oscillate around the MS on timescales of 
~ 0.4 tn, where tn is the Hubble time at that epoch. This 
corresponds to ~ 0.6 — 1.3 Gyr at 2 = 4 — 2. The more 
evolved galaxies, of higher masses at a given redshift, show 
fewer oscillations, one or two major compaction events fol¬ 
lowed by a decisive quenching. We investigate this further by 
performing a Fourier analysis in Section 5.5. Furthermore, 
we explore in Appendix B how the results based on a higher 
temporal resolution of snapshots compare with the results 
based on our standard resolution in Appendix B. In Fig¬ 
ure Bl, we plot the evolution of six simulated galaxies (11, 
12, 14, 25, 26, and 27) along the MS with high (upper pan¬ 
els) and standard (lower panels) temporal resolution. This 
figure confirms that the standard temporal resolution (steps 
of 100 Myr) traces the main features as well as the short 
term fluctuations. 

All galaxies evolve through certain sub-MS phases. At 
high redshifts, when the galaxy lives in a halo of relatively 
low mass (Mvir < 10 ^ Mq), these are only quenching at¬ 
tempts, i.e. the galaxy depletes some of its gas content and 
remains below the MS ridge for some time, but it quickly 
regains its position near or even above the MS ridge. This 
can be explained by fresh gas that flows through the halo 
into the galaxy, replenishing the gas in the disc, leading 
to a new compaction event and renewed star formation. 
Only after the galaxy’s halo reaches a critical virial mass of 
Mvir ^ 10 ^^'® Mq, and at late enough redshifts, the quench¬ 
ing process may continue for many 100 Myr up to several 
Gyr, allowing a successful quenching as opposed to a short¬ 
term quenching attempt at high- 2 . 

Focusing first on the colour gradients in Figure 4, we 
find that all eight galaxies show a strong positive correlation 
between Ams and Pgas.i. Galaxies on the upper envelope of 
the MS have about an order of magnitude higher central 
densities than galaxies below the MS, i.e., they are star¬ 
forming, compact blue nuggets, in which the central stellar 
density is soon to reach its peak value. Figure 5 shows that 
galaxies on the upper envelope of the MS tend to have a 
higher gas to stellar mass ratio than galaxies at the lower 
envelope. In Figure 6 , we see that galaxies at the top of the 
MS typically have a short tdep of ~ 300 Myr. On the other 
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Figure 4. Evolution of galaxy properties along and across the MS for eight simulated galaxies. Each panel shows the distance from the 
MS Ams = logio (sSFR/sSFRms) a-s a function of M*. The number in the upper left corner shows the identification number of the galaxy 
(see Table 1). The upper panels show four examples of low-mass galaxies, while the lower panels show examples of high-mass galaxies. 
The stars mark the redshifts z = 6, 4, and 2 from left to right. The colour coding corresponds to the gas density within the central 1 kpc, 
Pgas,l I corrected for the systematic dependence on z and M* using f{z) = § -I- C ^ + A and g{Mi,) = 77 -|- 7 X (log]^Q(M*) — 10.5), 

see Section 5.2 and Appendix C. The grey vertical line indicates the stellar mass at the time when the galaxy’s halo mass has reached 
Myij = 3 X 10^^ Af©, above which quenching is expected to be more likely. Galaxies at the top of the MS have about an order of 
magnitude higher Pgas.i than at the bottom of the MS. 


hand, galaxies at the bottom of the MS tend to have a longer 
tdep of np to ~ 1 Gyr. Qnenching galaxies all have low values 
of Pgas,i (< 10 ® M 0 /kpc®), low valnes of /gs (< 0 . 1 ), and 
long depletion times tdep (> 2 Gyr). 


5.2 Overall Population 

5.2.1 Gas-related gradients 

The trends in Figures 4, 5 and 6 are studied more quantita¬ 
tively in this section. We extend the focus from a few galax¬ 
ies to all the 26 simulated galaxies of this study. Figure 7 
shows from top to bottom the central gas density within 1 
kpc (pgas,i), the total gas mass (Mgas), gas to stellar mass 
ratio (/gs), and depletion time (tdep), corrected for the 2 - 
dependence and M*-evolution, as a function of the distance 
from the MS ridge (Ams) for all galaxies at 2 = 1 — 6 . 

We have corrected for the systematic time evolution of 
galaxies near the MS ridge because these key quantities may 
also evolve with cosmic time. In the simulations, we fol¬ 
low individual galaxies through cosmic time. This leads to 
a degeneracy between the redshift and the mass evolution, 
since all galaxies increase their mass with decreasing red- 


shift. To correct for the cosmic time evolution, we first fit 
the 2 -evolution of the galaxies that are near the MS ridge, 
I Ams I < 0.15. Varying the threshold for |Ams| between 
0.05 and 0.20 dex has no noticeable effect on the fits. Af¬ 
ter correcting for the 2 -evolution, we fit the M*-dependence. 
In this procedure, we assume that the cross-terms between 
2 -evolution and M*-dependence are negligible, i.e., we can 
practically determine them independently of each other. 

In Appendix C, Figure Cl, we show a more extended 
version of Figure 7. In the left panels, we show the uncor¬ 
rected quantities. In the middle-left and middle-right panels, 
we show the 2 -dependence and M*-dependence, respectively. 
The most right panels show the corrected quantities, i.e., 
the same as shown in Figure 7. We find a steep redshift evo¬ 
lution (/(z) = ^ -b ^ X logiQ(l -b 2 )) for Mgaa and /g^ with 
^ = —171±0.14 and -bl.63±0.17, respectively. For Pgaa.i and 
tdep, the 2 -evolution is much shallower with = -b0.30±0.26 
and —0.39 ± 0.19, respectively. Table 2 lists the best-fitting 
values. Following the 2 -correction, we fit the M*-dependence 
with the relation g(M*) = Ty-byx (log j^q(M*) —10.5). We find 
7 = +0.36±0.05, -b0.32±0.03, -0.15±0.03, and -0.19±0.03 
for pg 

as,l 5 fUgas, flfgas/AI©, and tciep, respectively (Table 2). 

In Figure 7, we fit each of the corrected quantities Q 
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Figure 5. Same as Figure 4, but the colour coding indicates the gas to stellar mass ratio /gg corrected to take out the overall 2 :-evolution 
and M*-dependence {f{z) and g{M^)). Galaxies at the top of MS tend to be more gas-rich, whereas galaxies at the bottom are typically 
gas-poor. 



2 -dependence 

M„- 

-dependence 

AMS"dependence 

^Pearson 


fR) = ^ + Cx logio(l-l-z) 

g{M„) = )7 ± 7 X (logjo M* - 10.5) 

Q = a + 5 X Ams 


This work 

C 

c 

V 

7 

a 

<5 


Pgas,l 

-1-7.28 ±0.12 

±0.30 ± 0.21 

±0.31 ± 0.05 

±0.36 ± 0.04 

-0.06 ±0.02 

±0.83 ±0.04 

0.57 ± 0.02 

+fgas 

±10.21 ± 0.09 

-1.71 ± 0.15 

±0.28 ±0.03 

±0.32 ± 0.03 

±0.00 ±0.01 

±0.48 ± 0.03 

0.53 ±0.02 

fgs 

-1.31 ± 0.06 

±1.63 ±0.10 

-0.13 ±0.02 

-0.15 ± 0.02 

±0.02 ±0.01 

±0.54 ±0.02 

0.70 ± 0.02 

^dep 

-0.19 ±0.07 

-0.39 ±0.11 

-0.17 ± 0.02 

-0.19 ±0.02 

±0.02 ±0.01 

-0.43 ± 0.02 

-0.51 ± 0.03 

Sm*,i 

±10.26 ±0.13 

-2.73 ±0.22 

±0.43 ±0.04 

±0.49 ± 0.04 

-0.06 ±0.01 

-0.12 ± 0.04 

-0.12 ± 0.02 

Re 

±0.70 ± 0.06 

-1.03 ±0.10 

-0.03 ±0.02 

-0.04 ±0.02 

±0.04 ±0.01 

±0.03 ±0.02 

0.05 ± 0.02 

n 

±0.75 ± 0.09 

-0.57 ± 0.14 

±0.12 ±0.03 

-0.12 ± 0.03 

±0.17 ± 0.01 

±0.16 ±0.03 

0.17 ± 0.02 

Observations" 

fgs 


±2.68 ± 0.05 


-0.37 ±0.04 


±0.49 ± 0.03 


^dep 


-0.34 ± 0.05 


±0.01 ± 0.03 


-0.49 ±0.02 



Table 2. List of best-fitting parameters for the gradients across the MS (Figures 7 and 8) in comparison with the values of observations. 
Last column shows the Pearson’s correlation coefficient rpearson- 

Data taken from Table 3 and 4 (global fits) of Genzel et al. (2015). 


with Q = a 8 x Ams? which is shown as a green line. We 
find 8 = +0.83 ± 0.03, 8 = +0.48 ± 0.02, 8 = +0.54 ± 0.02, 
and S = —0.43 ± 0.02 for Pgas.i, Mgas, fgs, and fdep, respec¬ 
tively (Table 2). Interestingly, the gradient of Pgas.i across 
the MS is the steepest, and the correlation is strongly posi¬ 
tive (Pearson’s correlation coefficient r = 0.57). This demon¬ 


strates that galaxies above the MS ridge have a significantly 
higher central gas density than galaxies below it: going from 
Ams = 0.5 dex below the MS ridge to 0.5 dex above it, we 
hnd that Pgas.i increases by nearly an order of magnitude. 
The gradients across the MS for the other quantities are very 
similar to each other (a slope of ~ 0.5) with |r| ~ 0.5 — 0.7. 
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Figure 6. Same as Figure 4, but the symbol size indicates the depletion time (tuep = A4gas/SFR) corrected to take out the overall 
2 -evolution and M*-dependence {f{z) and g{Mi,)). Galaxies at below the MS ridge tend to have longer depletion times than galaxies 
above it. 


Putting it all together, as seen by comparing Figure 4 to 
Figures 5 and 6 , the gradients across the MS are intimately 
related to the galaxy evolutionary phases of compaction to a 
blue nugget followed by central gas depletion and quenching. 
Galaxies above the MS are compact, gas-rich systems with 
high SFRs and short depletion times, whereas galaxies below 
are either more diffuse or depleted from central gas, of lower 
SFR and longer depletion times. 

5.2.2 Stellar-structure gradients 

In this section, we investigate stellar-structure gradients 
across the MS. We follow the approach used before for the 
gas-related gradients. Figure 8 shows from top to bottom the 
central stellar mass surface density within 1 kpc (Em*,i), the 
half-mass radius (Re), and the Sersic index (n), corrected 
for their systematic variations with redshift and with stellar 
mass, as a function of distance from the MS ridge (Ams) 
for all galaxies and snapshots at 2 = 6 — 1. As before, we 
show in Appendix C, Figure C2, a more extended version 
of Figure 8 , including the uncorrected quantities, and how 
we gradually account for the 2 -evolution and then the M*- 
dependence. The best-fitting values for the dependences on 
2 and on M*, as well as the gradients across the MS are 
listed in Table 2. 

In Figure 8 , the green lines show the linear regressions 
of the corrected quantities. They are much flatter in com¬ 
parison with the gas-related gradients from above. We find 


5 = -0.12±0.04, 5 = -b0.03±0.02, and <5 = -b0.16±0.03 for 
Svf*,!, Rb, and n, respectively. Furthermore, all three stellar 
structure quantities are not significantly correlated with the 
distance from the MS: the Pearson’s correlation coefficient 
is |r| < 0 . 20 . 


5.3 Gradients across the MS in Observations 

Genzel et al. (2015, hereafter G15) present CO-based and 
Herschel dust-based scaling relations of tdep and of /ga as a 
function of redshift, Ams and M*, for each of ~ 500 SFGs 
between 2 ~ 0 and 3. The best fit relations are spelled out 
in Table 2. 

Focusing first on /ga, we find a very similar scaling 
with 2 and M* in the simulations: G15 find as best fit 
for the 2 -dependence, combining CO and dust data^, = 
-1-2.68 ± 0.05 (-1-2.71 for the CO data and -1-2.32 for the 
dust data, respectively), which is slightly larger than our 
value of = -|-1.63±0.10. For the M*-dependence, G15 find 
7 = —0.37 ± 0.04 (—0.35 for the CO data and —0.40 for the 
dust data, respectively), also slightly larger than our value 


^ For the global combined CO-|-dust fit, G15 first added 0.1 dex 
to all GO values, and likewise subtracted 0.1 dex for all dust 
values before carrying out the global fit, in order to bring the two 
data sets to the same zero point. The values which we quote in 
the brackets are for the individual fits to the GO and dust data. 
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Figure 7. Gas-related galaxy properties across the MS. Shown 
from top to bottom are the gas mass within 1 kpc, total gas mass, 
gas mass to stellar mass ratio, and depletion time as a function 
of the distance from the MS Ams- The quantities are corrected 
for the systematic dependence on z and M*. Linear regression 
lines are shown (green lines), and the correlation coefficients are 
quoted. The color coding corresponds to stellar mass. The red 
dashed line and the shaded region show the best-fitting and its 
uncertainty of the observations by Genzel et al. (2015), indicating 
an excellent agreement with the gradients in the simulations. For 
the uncorrected quantities as well as the systematic dependence 
on redshift and stellar mass see Figure Cl. 
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Figure 8. Stellar-structure properties across the MS. Shown from 
top to bottom are the surface mass density within 1 kpc (Sjvf^ i), 
the half-mass radius (Ro), and the Sersic index (n), corrected for 
their systematic dependence on redshift and stellar mass, ver¬ 
sus the deviation from the MS ridge (Ams)- Linear regression 
lines are shown (green lines), and the correlation coefficients are 
quoted. The color coding corresponds to stellar mass. We see no 
significant systematic gradients across the MS. For the uncor¬ 
rected quantities as well as the systematic dependence on redshift 
and stellar mass see Figure C2. 


of 7 = —0.15 ± 0.02. However, despite these different scal¬ 
ings with z and M*, the obtained gradient across the MS 
is in very good agreement: G15 determined the gradient of 
/ga to be 5 = 0.49 ± 0.03 (0.53 for the CO data and 0.36 
for the dust data, respectively), while we measure a value of 
5 = 0.54 ±0.02. In Figure 7, the best-fitting line of G15 and 
its uncertainty are indicated as a red dashed line and as a 
red shaded area. 

Focusing now on tdep, we find an even better agree¬ 
ment between the ^-dependence and M*-dependence of G15 
and ours. For the ^-dependence, G15 find ^ = —0.34 ± 0.05 
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(—0.20 CO, —0.77 dust), while we find = —0.39 ± 0.11. 
For the M*-dependence, G15 find 7 = +0.01 ± 0.03 (—0.01, 
0.00), while we find 7 = —0.19 ± 0.02. For the tdep gradient 
across the MS, G15 determines S = —0.49 ± 0.02 (—0.43, 
—0.59), while we find S = —0.43 ± 0.02, i.e., our values are 
consistent within the systematic uncertainty. 

Wuyts et al. (2011) analyse the dependence of galaxy 
structure (size aud Sersic index) on the position of the galax¬ 
ies with respect to the MS ridge at 2 ; = 0 — 2. They per¬ 
formed the structural measurements on the light at the 
longest wavelength, high-resolution imaging, i.e., on zgso, 
7814 , and Tfieo, which are rest-frame UV to optical at z ~ 2. 
At 2 ~ 2, they find no significant gradient of galaxy size 
across the MS (Ams = —1.0 to 1.0). When galaxies lie more 
than two orders of magnitudes below the MS, i.e., quies¬ 
cent galaxies, the sizes decrease by ~ 0.2 dex. The Sersic 
index also tends to be roughly the same, n ~ 1 , between 
Ams = — 1.0 and 1 . 0 , i.e., there is no significant gradient of 
n across the MS. Again, when the galaxies are two orders 
of magnitude below the MS ridge, the Sersic index increases 
to n « 4.0. Our simulatious show the same null gradients in 
size and Sersic index about the MS ridge. As shown in Tac¬ 
chella et al. (2015a), for simulated galaxies that evolve along 
the MS ridge, the Sersic index increases from n ~ 1 — 2 at 
early times and low stellar masses (z < 2.5, M* < 10^° Mq) 
to n ~ 4 at later times and higher masses. Massive galax¬ 
ies that leave the MS have typically a high Sersic index of 
n ~ 4. In comparison to observations, where also the most 
massive galaxies on the MS have n ~ 2 — 3, it is important 
to consider that these measurements were performed on the 
UV and optical light, i.e., not on the mass as done in the 
simulations. Light-based profiles are indeed shallower, i.e., 
have a lower Sersic index, as the mass-based profiles, largely 
due to star-forming clumps in in the outskirts (Garollo et al. 
2014; Tacchella et al. 2015b). Furthermore, since we do not 
trace many galaxies well below the MS, we cannot address 
here totally quenched galaxies. 

Overall, we find excellent agreement between observa¬ 
tions and our simulations for the gradients across the MS, 
despite the different gas fraction estimates for a given M* 
(as discussed in Section 3). This indicates that the overall 
trends across the MS are robust, and that the intimate con¬ 
nection with the evolution through compaction, depletion, 
and quenching, and through phases of blue nuggets, is qual¬ 
itatively solid. 


5.4 Driver of the MS Gradients 

To understand the aforementioned gradients across the MS, 
we study the gas flow in the galaxies in more detail. In 
particular, we focus on the balance between the gas inflow 
rate (input term) and SFR plus gas outflow rate (drainage 
term). We expect galaxies that move upward on the MS, 
i.e., towards the blue nugget phase, are in an episode of 
compaction, where the central gas mass increases due to a 
high inflow rate and low outflow rate. On the other hand, in 
the post-compaction phase, we expect the opposite, namely, 
that the inflow rate is outweighed by SFR plus outflow rate. 

We investigate in Figure 9 the relation between the rate 
of change of distance from the MS, Ams, and the balance 
between the input term and drainage terms of gas in the cen- 



I 0 SIO ^5kpc 


logio 


inflow rate 
outflow rate + SFR 


Figure 9. In the main, large panel, we plot the rate of change 
of the distance from the MS, AmSi as a function of the bal¬ 
ance of gas input and drainage within 5 kpc, logiQ(B 5 i 5 p(.) = 
logio(inflow rate/(SFR + outflow rate)). The colour coding cor¬ 
responds to the position on the MS, Ams- We see a significant 
correlation between Ams and logpQ Bskpc with rpearson = 0.54. 
The solid line shows the best fit with Ams = (0.00 ± 0.07) + 
(5.5±0.3) X logjQ(iJ 5 kpc). This shows that galaxies that are mov¬ 
ing up towards the upper edge of the MS are inflow-dominated 
(compaction phase) whereas galaxies that are moving down to¬ 
wards the lower edge of the MS are depleted due to star formation 
and outflows while the inflow is suppressed (central quenching 
phase). In Figure D1 in Appendix D we split the balance term 
logic(T^skpc) inflow rate, outflow rate, and SFR. 


tral 5 kpc, namely logjQ Rskpc = logjQ[inflow rate/(SFR + 
outflow rate)]. Varying the radius considered between 3 and 
10 kpc does not change the result significantly. We have cho¬ 
sen 5 kpc as our fiducial scale since it best captures what 
happens within the galaxies towards their centres. A too 
large radius would only capture the gas exchange between 
the halo and the galaxy, and a too small radius would not 
capture whole extend of the inflow towards the galaxies’ 
centres. In Tacchella et al. (2015a), where we focus on the 
evolution of the surface density profiles, we find that the gas 
cusp of the compaction phase can reach out to 2 — 3 kpc. 
We have therefore chosen a slightly larger radius. Figure D1 
in Appendix D we split the balance term logjQ(i 35 kpc) into 
inflow rate, outflow rate, and SFR, each measured within 
the central 5 kpc. 

During ~ 70% of the time, this quantity is close to 0 
(i logic ^ 5 kpc| < 0.3), i.e., inflow rate « SFR + outfiow rate. 
However, there are episodes (lasting about 20% of the time) 
which are inflow-dominated (log^g iJskpc ^ 0.3 dex). The 
central gas density Pgas,i increases quickly (< 300 Myr) 
from 10^ Mq to a few times 10® Mq, i.e., this corresponds 
to dissipative, quick compaction phases of the galaxy gas 
into a compact, star-forming blue nugget (Zolotov et al. 
2015). Dekel & Burkert (2014) addressed the formation of 
blue nuggets by wet compaction using as an example the 
contraction associated with VDI. They applied the require¬ 
ment that for the inflow to be dissipative and therefore in- 
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tense, the characteristic timescale for star formation should 
be longer than the timescale for inflow, namely correspond¬ 
ing to logjg -Bskpc > 0. Otherwise, most of the disc mass will 
turn into stars before it reaches the bulge, the inflow rate 
will be suppressed, and the galaxy will become an extended 
stellar system. 

In the inflow-dominated phase, when the inflow rate 
is outweighing the SFR plus outflow rate in the centre, the 
rate of change of the distance from the MS, Ams, is positive, 
i.e., galaxies are moving up in the Ams — Af* plane of the 
universal MS (Figure 2). Most of these galaxies are below 
or just on the MS (Ams ^ 0). There are a few snapshots 
(< 3%) which are inflow-dominated, but where Ams < 0. 
In these cases, the galaxies are on the upper envelope of the 
MS (Ams ~ 0.3), i.e., they have reached the peak in the 
universal MS plane and are starting to move downwards. 

The inflow-dominated phases are followed by phases 
where log^g Bskpc ^ —0.3 dex, where the high SFR and 
the strong outflows, driven by the high SFR and stellar 
feedback, outweigh the inflow. This sudden suppression of 
the inflow at the end of the compaction process happens 
when the gas disc has shrunk and has not yet been replen¬ 
ished. This is the onset of a central depletion and therefore 
quenching phase, where the galaxies fall below the MS ridge 
(Ams < 0). As mentioned before, in low mass haloes at suf- 
hciently high redshifts, these are only quenching attempts, 
since gas quickly replenishes the disc. This gas is then avail¬ 
able to be triggered into a new episode of compaction and 
high SFR, which causes a subsequent quenching event, and 
so on. Full quenching of up to several Gyr into low sSFR 
signihcantly below the MS can be achieved preferentially at 
late redshifts, when the replenishment time is longer than 
the depletion time, and in particular after the galaxy’s halo 
reaches a critical virial mass of Mvir ^ Mq, corre¬ 

sponding to a stellar mass of ~ 10®'® Mq.. 

Overall, we find that the gas input and drainage within 
5 kpc is strongly correlated with the rate of change of the 
distance from the MS (rpearson = 0.54). From Figure D1 in 
Appendix D, we see that the individual rates are correlated 
to a lesser degree than the combined balance term Bskpc- 
The outflow rate and SFR are moderately correlated with 
Ams (rpearson = —0.32), while the inflow rate is not corre¬ 
lated with Ams- 


5.5 Oscillation Timescale 

As discussed before, galaxies oscillate about the MS ridge. 
To estimate the characteristic timescales for these oscilla¬ 
tions more quantitatively, we define the Fourier transform 
of the displacement from the MS ridge for each of the galax¬ 
ies individually: 

Aams, = ^ ^ ^MS„ exp f-27ri^) , (6) 

where n is the total number of snapshots available for a given 
galaxy in the redshift range z = 7 — 3. 

The median of the resulting Fourier spectrum over all 
galaxies is shown in Figure 10 as a black solid line. Also 
shown as dashed-blue and orange lines the are low- and high- 
mass galaxies, respectively. We find that the main period is 


t/tn 

1.55 0.55 0.30 0.19 0.14 0.11 



Figure 10. Fourier spectrum of the time for oscillations about 
the MS. The black line shows the median Fourier spectrum of 
all simulated galaxies and the shaded area marks the 16 and 84 
percentiles. The blue and orange dashed lines refer to the low- 
mass and high-mass subsamples. The green dashed line refers 
to the six galaxies for which we have thin time-steps available. 
The x-axis at the bottom is in units of inverse scale factor a~^, 
whereas the top axis is in units of cosmic time. The broad, global 
peak indicates that the timescale for the dominant oscillation is 
~ 0.2 — 0.5 til. The peak for the low-mass galaxies at 0.43 tp is 
narrower. 


~ 0.2 —0.5 tu for a full oscillation for galaxies evolving along 
the MS. The low-mass galaxies show a clear, narrow peak at 
~ 0.43 tu- This corresponds to ~ 0.4 Gyr at 2 = 6, and to 
~ 0.9 Gyr at 2 = 3. The spectrum for the high-mass galaxies 
is broader with a robust peak in the range 0.2 — 0.5 tu- 
The green dashed line in Figure 10 shows the Fourier 
spectrum for the high temporal resolution. We see that there 
is more power on small timescale fluctuations, as expected, 
but the peak at ~ 0.2 — 0.5 tu is confirmed. This is much 
longer than the dynamical time of the galaxy which is of 
the order of 30 Myr, i.e., the star-formation recipe or the 
feedback prescription used in these simulations do not drive 
the evolution about the MS ridge. 


6 THE ORIGIN OF CONFINEMENT OF THE 
MAIN SEQUENCE 

In this section, we explore the mechanisms responsible for 
confining the MS to a narrow strip about the MS ridge, con¬ 
necting the characteristic evolution pattern of high -2 galax¬ 
ies with the evolution along the MS. We apply a toy-model 
understanding for estimating the timescales that are en¬ 
coded into the MS. 
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6.1 Galaxy Properties across the universal MS 

We have seen in Sections 5.1 and 5.2 that differences in the 
position about the MS (i.e. in sSFR) at constant 2 ; and M* 
are associated with variations both in gas to stellar mass ra¬ 
tio (/gs oc and in depletion timescale (tdep oc 

in agreement with observational works by Magdis et al. 
(2012), Sargent et al. (2014), Huang & Kauffmann (2014), 
Silverman et al. (2015), Scoville et al. (2015), and G15. 
Galaxies above the MS ridge have larger gas fractions and 
smaller depletion times than galaxies at or below the MS 
ridge. 

The position about the MS has an even stronger depen¬ 
dence on the central gas density (pgas.i oc A^l^), i.e., on how 
the gas is distributed within the galaxies. This may indicate 
that the central gas density is the key factor involved in de¬ 
termining the MS width, i.e., an internal property has an 
important role in the confinement of the MS (though it may 
be stimulated externally, e.g. by a merger). The high central 
gas densities in more compact SFGs at the top of the MS 
lead to a decrease in the free-fall time tg, which itself leads 
to a shorter depletion time (tdep ~ tg/ea) even for a fixed 
total gas mass. Already Elbaz et al. (2011), Wuyts et al. 
(2011), Lada et al. (2012), and Sargent et al. (2014) put for¬ 
ward the idea that the decrease in tdep above the MS may 
be associated with internal parameters such as the central 
gas density. Furthermore, such compact SFGs (i.e., the blue 
nuggets) have been observationally detected (Barro et al. 
2013, 2014; Nelson et al. 2014; Bruce et al. 2014; Williams 
et al. 2014, 2015). 


6.2 Turnaround at the Upper and Lower Edge of 
the MS 

As highlighted in Section 5, galaxies oscillate around the MS 
equilibrium on timescales of 0.4 tn- Figure 11 sketches the 
evolution of a typical MS galaxy that eventually quenches 
its star formation at a late time in the massive end. During 
the oscillating evolution along the MS, a galaxy reaches min¬ 
ima and maxima in the distance from the MS (Ams). The 
massive galaxies (as measured at a given time, say 2 = 2) 
typically have one maximum, while less massive galaxies can 
have more than one. 

As shown before, the climb of a galaxy towards the top 
of the MS is due to a wet gas compaction, during which the 
gas inflow rate to the centre is faster than the SFR (Dekel 
& Burkert 2014). This leads to high SFR in a compact SFG 
(blue nugget), with high gas fraction and short depletion 
time. The turnaround at the top of the MS is a natural re¬ 
sult of more efficient gas depletion (shorter depletion time) 
by the high SFR and the associated feedback-driven out¬ 
flows, as observed by Gicone et al. (2016), combined with 
the suppression of gas inflow into the centre because the 
gas disc has shrunk and (at least temporarily) disappeared, 
or became gravitationally stable (morphological quenching, 
Martig et al. 2009; Genzel et al. 2014; Tacchella et al. 2015b). 

However, not all galaxies fully deplete and quench af¬ 
ter the first turnaround at the top of the MS. In galaxies of 
relatively low stellar mass, the low-mass halo allows rapid re¬ 
plenishment of the disc by fresh gas, with the replenishment 
time being shorter than the depletion time (trep < tdep). 


This sets the condition for another wet gas compaction, 
which can be triggered by mergers, counter-rotating streams 
or recycled gas, and can be associated with VDIs. The galax¬ 
ies, now at the lower envelope of the MS after the quenching 
attempt, turn around towards a new compact blue nugget 
phase with high SFR at the top of the MS. Once the re¬ 
plenishment becomes inefficient compared to the depletion 
(trep > tdep), typically at late redshifts and when the halo 
mass is above the threshold for virial shock heating, the con¬ 
ditions for wet compaction are not recovered. This allows the 
galaxy to quench inside-out all the way and thus drop below 
the MS. 

Self-regulation is the emerging feature that explains the 
small scatter about the MS ridge. Galaxies are not able to 
shoot above the MS ridge more than a few tenths of dex 
because the blue nugget phase naturally triggers central de¬ 
pletion, as the gas supply from the disc has been suppressed 
and the central gas is rapidly consumed by SFR and out¬ 
flows. On the other hand, at the lower envelope of the MS, 
the gas inflow quickly recovers, especially at high redshifts 
and in low halo masses, giving rise to a new compaction 
episode and an increase in star formation. Furthermore, the 
depletion and quenching at late cosmic times (2 < 2) and 
above a critical mass, with no push-back upwards through 
compaction, explains the downward slope of the MS at the 
high-mass end (e.g., Elbaz et al. 2007; Whitaker et al. 2014; 
Schreiber et al. 2015). 

6.3 Timescales for the Evolution on the MS 

The goal of this section is to estimate the timescales that 
are key to understand the evolution of galaxies along the MS 
and the confinement to it, namely the depletion time (tdep) 
and the replenishment time (trep). 

6.3.1 Depletion Time 

As discussed in Section 5.2 and shown in Figure Cl, in the 
simulations near the ridge of the MS, we measure an average 
depletion time of 

tdep = 0.44x(l + 2 )-“-x(^^^j^j Gyr. (7) 

This value and its dependence on cosmic time are consistent 
with observational estimates at 2 = 1 — 3 (Tacconi et al. 
2010, 2013; Genzel et al. 2015). 

The weak dependence of tdep on redshift, which is rather 
surprising at a first glance given the general growth of galac¬ 
tic dynamical times as (1 4 - 2 )“^'^^, may be qualitatively un¬ 
derstood by the compaction events as follows. As mentioned 
in Section 3, with a constant SFR efficiency eg, the varia¬ 
tion in tdep mostly reflects variations in tg (Krumholz et al. 
2012). In the Toomre regime, valid at high redshift, star for¬ 
mation occurs mostly in giant clumps. In these clumps, for 
a constant spin parameter for haloes and galaxies (in mass 
and redshift), one expects a systematic growth of tgep with 
cosmic time: 

tdep OC tg OC td OC tn OC (1 -b (8) 
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Compaction: 

Triggered by an intense gas infiow event, 
invoiving minor mergers or counter-rotating 
streams, and is commoniy associated with 
vioient disc instabiiity. The infiow rate is 
more efficient than the SFR. 
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Associated with a compact, massive core of gas and 
star-formation rate, short depietion time and high 
gas fraction. The downturn at the upper bound is 
due to the peak in SFR and outfiow and the 
suppression of infiow. Onset of quenching inside-out 
due to centrai gas depietion. 



Quenching Attempt: 

Centrai gas depietion gives rise to inside- 
out quenching. In low-mass haloes at high 
redshift, when trep < tdep, inflow of gas 
resumes, leading to the prerequisite for 
another compaction. 
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Quenching: 

Inflow rate cannot 
recover in hot haloes 
and at late times (frep 
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depletion and full 
quenching. 
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Figure 11. Sketch of the self-regulated evolution along the MS. SFGs are confined to a narrow MS before they quench. During this 
evolution, the galaxy lives through one or more blue nugget phases during which a minimum in gas depletion time and a maximum in gas 
fraction are reached. The blue nugget phases are followed by gas depletion inside-out. This quenching attempts fail for low halo masses 
and at high redshifts since the recovered inflow rate triggers a new episode of compaction and high star formation. At high halo masses 
(hot halo), the inflow rate cannot recover and the galaxies ceases its star formation activity. 


where td is the galaxy dynamical crossing time. This is in 
contrast with the slow growth of tdep in the simulations and 
observations. 

At low redshifts, as argued in Krumholz et al. (2012), 
there is a transition to the giant molecular cloud regime, 
where the surface density is constant. Then, the growth of 
tff with time is suppressed, and one expects a similar effect 
on tdep- However, it turns out that the predicted slowdown 
in the growth rate of tdep occurs too late and is insufficient 
for explaining the indicated slow growth of tdep- 

Another possibility is that the sequence of wet com¬ 
paction events may provide a clue for the suppressed growth 
of tdep- Each such event causes a decrease in the dynami¬ 
cal time of the galaxy, and thus a corresponding decrease 
in tdep, balancing the natural systematic increase in time 
based on Equation 3. This can be addressed in conjunction 
with the oscillations in the MS- When passing through the 
ridge on the way up (increasing Ams) at a later time, at 
the early stages of a compaction process, the system is still 


extended, so tdep oc td is still relatively long, roughly fol¬ 
lowing Equation 8- However, when passing through the MS 
ridge on the way down (decreasing Ams), during the early 
stages of the post-compaction quenching process, the sys¬ 
tem is still gas-rich, with a short tdep oc td- The compaction 
causes a decrease in td, which balances the natural increase 
in time from Equation 8- At high-z, there are galaxies mov¬ 
ing both up and down the MS- However, at later «, more 
galaxies are at their post-blue-nugget phase, moving down 
at the MS ridge, allowing the compaction-driven decline of 
td to balance the natural cosmological growth of td- Since 
more massive galaxies quench earlier, and to higher densi¬ 
ties, we expect their tdep at the MS ridge to become shorter, 
and at earlier times, as seen in the simulations- 

Our attempt to address the origins of the time evolution 
of tdep is only a qualitative preliminary step- What matters 
for the arguments below concerning the confinement of the 
MS is the general evolution of tdep with cosmic time and 
mass. 
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redshift 


Figure 12. Sketch of expected quenching efficiency in the plane 
of halo mass Mvir and redshift 2 . Colour coding and contours 
correspond to the expected quenching efficiency, which is given 
by the ratio of the replenishment time trep (Equation 9) and the 
depletion time fdep (Equation 7). The horizontal line at Afvir = 
Mq crudely marks the threshold mass for a stable shock 
based on spherical infall analysis (Dekel &: Birnboim 2006), which 
in practice stretches over an order of magnitude in mass. Below 
this curve the flows are predicted to be predominantly cold and 
above it a shock-heated medium is expected to extend out to the 
halo virial radius. The inclined dashed curve is the conjectured 
upper limit for cold streams in hot haloes, valid at redshifts higher 
than -Zcrit ~ 2. 


6.3.2 Replenishment Time 

After the central depletion started at the blue nugget phase, 
it can either bounce back at the lower edge of the MS into a 
new compaction phase, or it can continue to quench to well 
below the MS (Figure 11). The critical times to be compared 
to the depletion time are the time for replenishment of the 
gas disc, and the time for the next intense accretion episode, 
e.g. a merger, that can trigger a new compaction event. As 
long as the halo is not massive and hot enough to suppress 
the streaming of cold gas through it, these two timescales 
can be approximated by the timescale for mass accretion into 
the galaxy, the inverse of the specific accretion rate given in 
Equation 4, namely, 

( l.GOGyr z = 2 

t,ep ~ 25so.j4(l + ® Gyr ~ i 0.78Gyr 2 = 3, (9) 

[o.45Gyr 2 = 4 

where so,o4 = Sh/0.04 « 1. 

The replenishment time can become much longer if 
the halo is more massive than a threshold mass, Mvir ~ 
Mshock ~ Mq, such that it can support a stable virial 

shock that keeps the circum-galactic medium at the virial 
temperature, and when cold streams are suppressed at late 
redshifts (Birnboim & Dekel 2003; Dekel & Birnboim 2006). 


6.3.3 Quenching Attempt versus Full Quenching 

As mentioned before, full quenching is achieved when the 
replenishment time is longer than the depletion time. In 
Figure 12, we show a very rough estimate of the quench¬ 
ing efficiency, dehned as the ratio of the estimated replen¬ 
ishment time trep (Equation 9) and the estimated depletion 
time tdep (Equation 7), in the plane of the halo mass Mvir 
and redshift z. The red line indicates the boundary where 
trep = tdep. We cautiou the reader that at 2 < 1, Equation 9 
gives only a very crude estimate of the replenishment time 

trep. 

While being above the red line where trep ~ tdep indi¬ 
cates that the quenching process can possibly proceed, the 
actual value of trep/^dep where full, long-term quenching is 
achieved may be larger, on the order of a few. For one thing, 
the newly accreted gas might not be available immediately 
for star formation, which can cause a delay. This delay might 
be larger for galaxies towards lower redshifts and it could 
also depend on stellar mass. Thus, in Figure 12, long-term 
quenching can be crudely expected in the red area above 

^rep/tdep 10. 

In the regime of Mvir < Mshock ~ 10^^’^, comparing 
Equation 7 with Equation 9, we see that, for a galaxy with 
a stellar mass of a few times Mq, the condition for 

full quenching is valid for 2 < 2.5. For massive galaxies, the 
condition is valid earlier (2 < 3.0), whereas for lower-mass 
galaxies, the condition is valid later (2 < 2.0). This may 
explain the decisive quenching of massive galaxies at high- 
2 , even before halo quenching dominates. In addition, the 
hot medium in haloes of Mvir > Mshock at 2 > 2crit ~ 2 is 
predicted to host penetrating cold streams, while haloes of a 
similar mass at 2 < 2crit are expected to be all hot, shutting 
off most of the gas supply to the inner galaxy. Therefore, 
once Mvir > Mshock, the hot halo can make trep much longer, 
so the condition trep > tdep is more easily fulfilled even at 
high redshifts (see also Fig. 7 of Dekel & Birnboim 2006). 

6.4 Further Estimates 

First, we estimate the timescale for compaction and check 
whether it is consistent with the Fourier analysis of the MS 
oscillations presented in Section 5.5. Secondly, we calculate 
the width of the MS based on several simple assumptions. 
We caution the reader that both the estimate for the com¬ 
paction timescale and the calculation of the width of the 
MS are rather crude and should be regarded as consistency 
checks. 


6 . 4.1 Compaction Time 

Galaxies do not remain “sub-” or “super-” MS galaxies, but 
they rather oscillate about the MS ridge. In the simulations, 
e.g.. Fig. 16 of Zolotov et al. (2015), the duration of the 
compaction phase from the onset of compaction to the blue 
nugget is on average about 

f l.OOGyr 2 = 2 

Aom ~ (0.3 - 0.4)tH =2 < 0.66Gyr 2 = 3, (10) 

[o.50Gyr 2 = 4 
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where tn is the Hubble time at the blue nugget phase. 

We can very crudely estimate the expected ballpark of 
tcom in the following two ways. If the compaction phase is 
driven by a minor merger event, we expect tcom to be compa¬ 
rable to the duration of a minor merger, from the first close 
passage to coalescence. This is in the ball park of the halo 
crossing time, tvir ~ 0.2 fn (based on spherical collapse), 
which is not far from what we see in the simulations. 

Alternatively, assuming VDI-driven wet inflow (Dekel 
et al. 2009a,b; Dekel & Burkert 2014), we can evaluate tcom 
by the evacuation time of the disc, 

tcom ~ ^ tmig ~ SOtd ~ 2tvir ~ O.dtH; (H) 

assuming a ~ 0.2 for the fraction of cold disc mass in 
clumps. The timescale tmig is the migration time of the 
clumps, and td is the disc crossing time, which can be es¬ 
timated by td ~ Atvir with the spin parameter A ~ 0.04. 
This estimate recovers again the ballpark of tcom from the 
simulations. 

6.4.2 Oscillation Timescale 

In Section 5.5, we measured the timescale for a full oscilla¬ 
tion at 2 = 3 — 6 based on a Fourier analysis to be about 
0.2 — 0.5 tn. We can now estimate this oscillation timescale 
from the derived timescales for compaction and replenish¬ 
ment. From Equations 9 and 10, we obtain for the oscillation 
timescale tosc ~ tcom + trep ~ (0.3 -I- 0.3) tu ~ 0.6 tn, which 
is roughly consistent with the Fourier analysis. 

6.4.3 The Width of the MS 

In Section 4, we determined ums = 0.27 at 2 = 3 — 6, which 
is in good agreement with observations. Based on the phys¬ 
ical picture above, we can attempt to estimate the expected 
scatter of the MS from the estimated timescales. 

Consider a galaxy that moves from the upper to the 
lower edge of the MS over trep ~ tdep. Assume that during 
this phase, the gas mass is about half its peak value at the 
blue nugget point on average, so the SFR is about half its 
peak value (by the Kennicutt law). Since there is not much 
inflow in this phase, the baryon mass is roughly constant. 
Therefore, if the gas mass is half what it was at the peak, 
the stellar mass is roughly twice what it was. This means 
that the sSFR has dropped by more than a factor of 4. Dur¬ 
ing the quenching episode, say between 2 = 4 and 3, the 
sSFR MS ridge has dropped by a factor of ~ 1.6. This gives 
CMS ^ logio(-\/4/1.6) ~ 0.20 dex, similar to the observed 
and simulated width of the MS. 

Next, consider a galaxy during its compaction phase, 
from the lower to the upper edge of the MS. During this 
phase, on average, the gas density in the centre increases by 
a factor of order 10, while the total gas mass is changing by a 
much smaller factor. Based on Equation 3, the SER increases 
by a factor of vTO ~ 3.2. The stellar mass may be doubling 
its value from the green to the blue nugget point. Therefore, 
the sSFR increases by a factor of ~ 1.6. As before, the sSFR 
MS ridge has dropped during the compaction phase by a 
factor of ~ 1.6. This gives cms > logio('\/l.6 x 1.6) « 0.20 


dex, again in the ball park of the observed and simulated 
width of the MS. 


7 DISCUSSION 

In this section, we discuss the implications of the picture 
of the MS as discussed above for the quenching of galaxies. 
Furthermore, we highlight how the outlined picture of the 
MS can be confirmed in observations. Finally, we highlight 
some caveats of the analysis presented here and how one can 
improve in future work. 


7.1 Implications for Quenching 

There is solid observational evidence that the cessation of 
star formation in some galaxies, which results in the emer¬ 
gence of quiescent galaxies, correlates with both galaxy mass 
and environment (e.g.. Dressier 1980; Balogh et al. 2004; 
Baldry et al. 2006; Kimm et al. 2009; Peng et al. 2010; Woo 
et al. 2013; Knobel et al. 2013; Kovac et al. 2014). Fur¬ 
thermore, quenching has been observed to correlate strongly 
with morphology and galaxy structure in the local universe 
(Kauffmann et al. 2003; Franx et al. 2008; Cibinel et al. 2013; 
Fang et al. 2013; Schawinski et al. 2014; Bluck et al. 2014; 
Woo et al. 2015) and at high -2 (Wuyts et al. 2011, 2012; Bell 
et al. 2012; Cheung et al. 2012; Szomoru et al. 2012; Barro 
et al. 2013; Lang et al. 2014; Tacchella et al. 2015c). On the 
other hand, quenched disc galaxies have also been observed 
(McGrath et al. 2008; van Dokkum et al. 2008; Bundy et al. 
2010; Salim et al. 2012; Bruce et al. 2012; Carollo et al. 

2014) . 

The physical nature of quenching, and its link to mor¬ 
phology, is a central issue in galaxy evolution. From theory, 
it has been proposed that galaxies starve out of gas by rapid 
gas consumption into stars, in combination with the associ¬ 
ated outflows driven by stellar feedback (e.g., Dekel & Silk 
1986; Murray et al. 2005) or super-massive black hole (e.g., 
Di Matteo et al. 2005; Croton et al. 2006; Ciotti & ©striker 
2007; Cattaneo et al. 2009), and/or by a slowdown of gas 
supply into the galaxy (e.g., Rees & ©striker 1977; Dekel & 
Birnboim 2006; Hearin & Watson 2013; Feldmann & Mayer 

2015) . Another possibility includes morphological quench¬ 
ing, which argues that the growth of a central mass concen¬ 
tration, i.e., a massive bulge, stabilizes a gas disc against 
fragmentation (Martig et al. 2009). 

Using observations alone, it is difficult to constrain 
the physical nature of quenching because a correlation be¬ 
tween quenching and a certain galaxy property does not 
necessarily imply a direct causal relation or the direction 
of such a causality (e.g., Carollo et al. 2014). For exam¬ 
ple the observed correlation between central stellar density 
and quenching could either arise because high stellar density 
causes quenching, or because another property that is asso¬ 
ciated with high stellar density causes quenching, or because 
quenching leads to high stellar density. The other property 
may be, for example, central gas density, AGN feedback, or 
halo mass. As an example, Lilly & Carollo (2016, in prepa¬ 
ration) show that a central surface density threshold for 
quenching could in principle be the result of a hierarchical 
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accretion of mass in galaxies, in which galaxies that quench 
earlier are denser. 

Tacchella et al. (2015b) mapped out the M* and SFR 
distribution on scales of 1 kpc in « ~ 2 SFGs. They fonnd 
that the ~ 10^^ Mq galaxies quench inside-out, where the 
star formation in the centre ceases within < 200 Myr after 
a ~ 2, whereas the outskirts still form stars for 1 — 3 Gyr. 
In our simulations, we find a similar inside-out quenching 
signature where the depletion starts from the centre in the 
post-blue-nugget phase (see also Tacchella et al. 2015a, for 
a detailed analysis of the evolution for the M* and SFR 
profiles in the simulations). As shown above, our simulated 
galaxies are able to deplete and quench rapidly at 2 ~ 1 — 3, 
forming galaxies that resemble today’s typical ~ 10^^ Mq 
galaxies, which show features of a gas-rich, dissipative for¬ 
mation process (e.g., Bender et al. 1988; Carollo et al. 1993; 
Faber et al. 1997; Cappellari et al. 2007). 

From the analysis of the simulations presented here, we 
can now make a step forward in understanding the phys¬ 
ical nature of quenching. Our simulations show that the 
interplay between the gas and the stellar components to¬ 
gether with the dark-matter halo is able to explain quench¬ 
ing of galaxies. Our analysis highlighted the key role of halo 
mass, which directly translates to the gas replenishment 
time of the galaxy’s disc (see Section 6.3.3). We have seen 
that galaxies depleted their centres first, i.e., quenching pro¬ 
gresses inside-out. These observational signatures have been 
found also in observations of 2 ~ 2 galaxies. Furthermore, 
the observation that the central stellar mass density (bulge 
mass) is a good indicator for quiescence (Kauffmann et al. 
2003; Franx et al. 2008; Cheung et al. 2012; Fang et al. 
2013; Bluck et al. 2014; Lang et al. 2014; Woo et al. 2015; 
Barro et al. 2015) is in good agreement with the thoughts 
presented here. If a compaction with a nuclear “starburst” 
precedes quenching, one expects a high central stellar mass 
density to have built up. Similarly, due to the key role of halo 
mass in our picture, and the correlation of stellar mass with 
halo mass, we would also expect that total stellar mass is a 
good measure for quiescence (e.g., Peng et al. 2010; Carollo 
et al. 2013). 


7.2 Observational Consequences 

As demonstrated above, our simulations reveal gradients of 
gas fraction and depletion time across the MS that are con¬ 
sistent with the observed gradients by G15, and explain the 
phenomena by the evolution through compaction, depletion, 
and quenching events. Our additional key prediction is a 
strong gradient of the core gas density, e.g. Pgaa.i. Future ob¬ 
servations that will resolve the gas distribution within indi- 
vidnal galaxies should be able to confirm this. Furthermore, 
star formation is expected to occur centrally concentrated 
at the top of the MS, whereas it is predicted to be ring-like 
distributed in massive systems at the lower envelope of the 
MS proceeding the blue-nugget phase before quenching. In¬ 
dications for this have already been seen in Genzel et al. 
(2014) and Tacchella et al. (2015b). 


7.3 Caveats 

We have mentioned several limitations of the simulations 
used here in Section 2.2. For example, our simulations do not 
include AGN feedback and feedback associated with cosmic 
rays and magnetic field. Nevertheless, Geverino et al. (2014) 
have shown that our simulations match basic observations, 
such as SFRs, gas fractions, and stellar-to-halo mass frac¬ 
tions, at least at least within observational uncertainties. If 
the SFR at very high redshifts is still overestimated, the dra¬ 
matic events in the evolution of galaxies that concern us here 
may occur somewhat earlier than in the real universe. How¬ 
ever, it seems that in the current simulations these evens, 
the compaction and the subsequent onset of quenching, do 
occnr at cosmological times that are consistent with obser¬ 
vations. 

The number of simulated galaxies analysed in this work 
is limited to 26, outputted in ~ 900 snapshots from 2 = 6—1. 
Glearly, a larger number of simulated galaxies would have 
been better for studying the average properties of galaxies 
as a function of time and mass, and the scatter about them, 
but the number of galaxies is limited being computationally 
expensive. However, the key aspect of this work is that we 
resolve the gas physics of high-density star-forming regions 
at ~ 25 pc scales in a cosmological context, which allows 
us to capture the evolution of galaxies along and about the 
MS. This is key to understand the evolution of galaxies on 
the MS. 


8 CONCLUSION 

Our zoom-in cosmological simulations of relatively massive 
galaxies at 2 > 1 reveal that the SFGs are confined to a uni¬ 
versal MS as they evolve through gas compaction, depletion, 
and replenishment. They reproduce the observed properties 
of the MS and allow us to present a simple understanding 
for the evolution and the scatter of the MS. 

In the simulations, in the redshift range 2 = 6—1, the 
MS ridge of sSFR declines in time, sSFRms oc (1-1-2)^'®, with 
only a weak mass dependence. This follows closely the pre¬ 
dicted decline of the specific halo-mass accretion rate in the 
Einstien-deSitter phase of the expanding universe. We find 
that galaxies oscillate about this MS ridge on timescales of 
~ 0.2—0.5 tn, which corresponds to 0.2—0.9 Gyr at 2 = 6—3. 
The SFGs are confined to a narrow MS of width ±0.27 dex. 
The simulated evolution with redshift, the MS width, and 
the bending of the MS at high masses are consistent with 
observations. 

The simulations reveal that the high-SFR galaxies at 
the upper envelope of the MS tend to be compact blue 
nuggets with high gas fractions and short depletion times. 
The lower-SFR galaxies at the lower envelope of the MS 
typically get there after partial central gas depletion that is 
naturally triggered at the blue-nugget phase, and they have 
lower gas fractions and longer depletion times. The mea¬ 
sured gradients of gas fraction and depletion times across the 
MS are in agreement with observations. The steepest gradi¬ 
ent across the MS is predicted for the central gas density, 
Pgas.i, which indicates that the main mechanism responsible 
for the width of the MS involves an internal property of the 
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galaxy, though it may be determined by an external process. 
On the other hand, we do not find any significant gradients 
across the MS for stellar structure properties, namely central 
stellar mass surface density, size, and Sersic index. 

Propagation up towards the upper envelope of the MS 
(blue nugget phase) is associated with gas compaction of the 
gas disc, triggered by a minor merger or counter-rotating 
streams through violent disc instability, where the inflow 
rate is higher than the SFR. The blue nugget phase marks 
the onset of gas depletion, during which the central gas is 
exhausted due to SFR and outflow, while the inflow from the 
disc that has shrunk is suppressed. In the post-compaction 
phase the galaxies show the signature of inside-out quench¬ 
ing, where the SFR first reduces in the centre. After reaching 
the bottom of the MS, an upturn can occur if the extended 
disc has been replenished by fresh gas and a new compaction 
event is triggered prior to total gas depletion, namely, if 
trep < tdep. The self-regulated nature of these mechanisms 
can explain the confinement of SFGs into a MS narrower 
than ±0.3 dex. 

Full quenching with a departure from the MS occurs in 
massive haloes or at low redshifts, where the disc replen¬ 
ishment is slow compared to the depletion time, which we 
find to be 0.2 — 1.0 Gyr, declining with mass and increasing 
with cosmic time, consistent with observations. The require¬ 
ment for full quenching, t^ep < frep, tends to be fulfilled at 
redshifts a Si 3, and possibly at higher redshifts for more 
massive galaxies. A halo above the critical mass for virial 
shock heating, ~ 10 ^^ ®Mq, helps suppressing the replenish¬ 
ment and causing long-term quenching. 


ACKNOWLEDGEMENTS 

We acknowledge stimulating discussions with Guillermo 
Barro, Sandy Faber, Reinhard Genzel, Mark Krumholz, Si¬ 
mon Lilly, Benny Trakhtenbrot, Joanna Woo, and Adi Zolo¬ 
tov. We thank the referee for constructive and useful com¬ 
ments, which helped to improve the manuscript. ST thanks 
AD and his group for the hospitality during his visit at 
HUJI. Development and most of the analysis have been per¬ 
formed in the astro cluster at HU. The simulations were 
performed at the National Energy Research Scientific Gom- 
puting Center (NERSC), Lawrence Berkeley National Lab¬ 
oratory, and at NASA Advanced Supercomputing (NAS) 
at NASA Ames Reserach Center. This work was supported 
by ISF grant 24/12, by GIF grant G-1052-104.7/2009, by 
the I-CORE Program of the PBG, ISF grant 1829/12, by 
MINECO grant AYA2012-32295, by GANDELS grant HST- 
GO-12060.12-A, and by NSF grants AST-1010033 and AST- 
1405962. We acknowledge support by the Swiss National 
Science Foundation. 

REFERENCES 

Abramson L. E., Gladders M. D., Dressier A., Oemler Jr. A., 
Poggianti B., Vulcani B., 2015, ApJ, 801, L12 [1] 

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, 
ApJ, 770, 25 [2.1] 


Baldry I. K., Balogh M. L., Bower R. G., Glazebrook K., Nichol 
R. C., Bamford S. P., Budavari T., 2006, MNRAS, 373, 469 

[7.1] 

Balogh M. L., Baldry I. K., Nichol R., Miller G., Bower R., Glaze- 
brook K., 2004, ApJ, 615, LlOl [7.1] 

Barro G., et ah, 2013, ApJ, 765, 104 [2.2, 6.1, 7.1] 

Barro G., et ah, 2014, ApJ, 791, 52 [6.1] 

Barro G., et ah, 2015, preprint, (arXiv: 1509.00469) [7.1] 
Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379 
[2.2, A] 

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57 [1, 
2.2, Al, A] 

Bell E. F., et ah, 2012, ApJ, 753, 167 [7.1] 

Bender R., Doebereiner S., Moellenhoff C., 1988, A&AS, 74, 385 

[7.1] 

Birnboim Y., Dekel A., 2003, MNRAS, 345, 349 [6.3.2] 

Birnboim Y., Dekel A., Neistein E., 2007, MNRAS, 380, 339 [1] 
Bluck A. F. L., Mendel J. T., Ellison S. L., Moreno J., Simard L., 
Patton D. R., Starkenburg E., 2014, MNRAS, 441, 599 [7.1] 
Bouche N., et ak, 2010, ApJ, 718, 1001 [1, 4.2] 

Bournaud F., Elmegreen B. G., Elmegreen D. M., 2007, ApJ, 670, 
237 [1] 

Bournaud F., et al., 2012, ApJ, 757, 81 [1] 

Brinchmann J., Chariot S., White S. D. M., Tremonti C., Kauff- 
mann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 
1151 [1] 

Bruce V. A., et ah, 2012, MNRAS, 427, 1666 [7.1] 

Bruce V. A., et ah, 2014, MNRAS, 444, 1660 [6.1] 

Bryan G. L., Norman M. L., 1998, ApJ, 495, 80 [2.3] 

Bundy K., et ah, 2010, ApJ, 719, 1969 [7.1] 

Burkert A., et ah, 2010, ApJ, 725, 2324 [1] 

Burkert A., et ah, 2015, preprint, (arXiv: 1510.03262) [2.2, Al, 

A] 

Cacciato M., Dekel A., Genel S., 2012, MNRAS, 421, 818 [1] 
Cappellari M., et ah, 2007, MNRAS, 379, 418 [7.1] 

Carollo C. M., Danziger I. J., Buson L., 1993, MNRAS, 265, 553 

[7.1] 

Carollo C. M., et ah, 2013, ApJ, 773, 112 [1, 7.1] 

Carollo C. M., et ah, 2014, preprint, (arXiv: 1402.1172) [5.3, 
7.1] 

Cattaneo A., et ah, 2009, Nature, 460, 213 [7.1] 

Ceverino D., Klypin A., 2009, ApJ, 695, 292 [2.1, 2.1] 

Ceverino D., Dekel A., Bournaud F., 2010, MNRAS, 404, 2151 

[ 2 . 1 , 2 . 2 ] 

Ceverino D., Dekel A., Mandelker N., Bournaud F., Burkert A., 
Genzel R., Primack J., 2012, MNRAS, 420, 3490 [2.1, 2.2] 
Ceverino D., Klypin A., Klimek E. S., Trujillo-Gomez S., 
Churchill C. W., Primack J., Dekel A., 2014, MNRAS, 442, 
1545 [2, 2.1, 2.2, 7.3] 

Ceverino D., Dekel A., Tweed D., Primack J., 2015a, MNRAS, 
447, 3291 [2.2] 

Ceverino D., Primack J., Dekel A., 2015b, MNRAS, 453, 408 [2] 
Chabrier G., 2003, PASP, 115, 763 [2.1] 

Cheung E., et ah, 2012, ApJ, 760, 131 [7.1] 

Cibinel A., et ah, 2013, ApJ, 776, 72 [7.1] 

Cicone C., Maiolino R., Marconi A., 2016, preprint, 

(arXiv; 1601.04715) [6.2] 

Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038 [7.1] 

Conroy C., Wechsler R. H., 2009, ApJ, 696, 620 [2.2, A] 

Croton D. J., et ah, 2006, MNRAS, 365, 11 [7.1] 

Daddi E., et ah, 2007, ApJ, 670, 156 [1, 4.2] 

Daddi E., et ah, 2010, ApJ, 713, 686 [1] 

Dave R., Oppenheimer B. D., Finlator K., 2011, MNRAS, 415, 

11 [ 1 ] 

Dave R., Finlator K., Oppenheimer B. D., 2012, MNRAS, 421, 
98 [1] 

Dayal P., Ferrara A., Dunlop J. S., 2013, MNRAS, 430, 2891 [1] 


MNRAS 000, 1-27 (2016) 


22 S. Tacchella, A. Dekel, C. M. Carollo, et al. 


Dekel A., Birnboim Y., 2006, MNRAS, 368, 2 [2.3, 5.1, 12, 6.3.2, 
6.3.3, 7.1] 

Dekel A., Burkert A., 2014, MNRAS, 438, 1870 [1, 5.4, 6.2, 6.4.1] 
Dekel A., Krumholz M. R., 2013, MNRAS, 432, 455 [2.2] 

Dekel A., Mandelker N., 2014, MNRAS, 444, 2071 [1, 4.2] 

Dekel A., Silk J., 1986, ApJ, 303, 39 [2.1, 7.1] 

Dekel A., et al., 2009a, Nature, 457, 451 [6.4.1] 

Dekel A., Sari R., Ceverino D., 2009b, ApJ, 703, 785 [1, 6.4.1] 
Dekel A., Zolotov A., Tweed D., Cacciato M., Ceverino D., Pri- 
mack J. R., 2013, MNRAS, 435, 999 [1, 4.1, 4.2] 

Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604 

[7.1] 

Dressier A., 1980, ApJ, 236, 351 [7.1] 

Dutton A. A., van den Bosch F. C., Dekel A., 2010, MNRAS, 405, 
1690 [1] 

Elbaz D., et al., 2007, A&A, 468, 33 [1, 4.2, 6.2] 

Elbaz D., et al., 2011, A&A, 533, A119 [1, 6.1] 

Faber S. M., et al., 1997, AJ, 114, 1771 [7.1] 

Fakhouri O., Ma C.-P., 2009, MNRAS, 394, 1825 [1] 

Fang J. J., Faber S. M., Koo D. C., Dekel A., 2013, ApJ, 776, 63 

[7.1] 

Feldmann R., 2015, MNRAS, 449, 3274 [1] 

Feldmann R., Mayer L., 2015, MNRAS, 446, 1939 [7.1] 

Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., King- 
don J. B., Verner E. M., 1998, PASP, 110, 761 [2.1] 

Forbes J. C., Krumholz M. R., Burkert A., Dekel A., 2014, MN¬ 
RAS, 443, 168 [1] 

Franx M., van Dokkum P. G., Schreiber N. M. F., Wuyts S., 
Labbe I., Toft S., 2008, ApJ, 688, 770 [7.1] 

Gammie C. F., 2001, ApJ, 553, 174 [1] 

Genel S., Bouche N., Naab T., Sternberg A., Genzel R., 2010, 
ApJ, 719, 229 [1] 

Genzel R., et al., 2010, MNRAS, 407, 2091 [1] 

Genzel R., et al., 2014, ApJ, 785, 75 [6.2, 7.2] 

Genzel R., et al., 2015, ApJ, 800, 20 [1, 5, 2, 5.3, 7, 6.3.1] 
Gladders M. D., Oemler A., Dressier A., Poggianti B., Vulcani 
B., Abramson L., 2013, ApJ, 770, 64 [1] 

Haardt F., Madau P., 1996, ApJ, 461, 20 [2.1] 

Hearin A. P., Watson D. F., 2013, MNRAS, 435, 1313 [7.1] 
Hopkins P. F., Keres D., Murray N., Quataert E., Hernquist L., 
2012, MNRAS, 427, 968 [2.2] 

Huang M.-L., Kauffmann G., 2014, MNRAS, 443, 1329 [1, 6.1] 
Kauffmann G., et al., 2003, MNRAS, 346, 1055 [7.1] 

Kelson D. D., 2014, preprint, (arXiv: 1406.5191) [1] 

Kennicutt Jr. R. C., 1998, ARA&A, 36, 189 [2.1] 

Kimm T., et al., 2009, MNRAS, 394, 1131 [7.1] 

Knobel C., et al., 2013, ApJ, 769, 24 [7.1] 

Komatsu E., et al., 2009, ApJS, 180, 330 [2.3] 

Kovac K., et al., 2014, MNRAS, 438, 717 [7.1] 

Kravtsov A. V., 2003, ApJ, 590, LI [2.1] 

Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 
73 [2.1] 

Krumholz M. R., Dekel A., 2010, MNRAS, 406, 112 [2.2] 
Krumholz M. R., Dekel A., 2012, ApJ, 753, 16 [2.2] 

Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155 [2.2] 
Krumholz M. R., Dekel A., McKee G. F., 2012, ApJ, 745, 69 [2.1, 
3, 6.3.1, 6.3.1] 

Lada C. J., Forbrich J., Lombardi M., Alves J. F., 2012, ApJ, 
745, 190 [6.1] 

Lang P., et al., 2014, ApJ, 788, 11 [7.1] 

Leitherer C., et al., 1999, ApJS, 123, 3 [2.1] 

Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, 
ApJ, 772, 119 [1, 3, 4.2] 

Magdis G. E., et al., 2012, ApJ, 760, 6 [1, 6.1] 

Magnelli B., et al., 2014, A&A, 561, A86 [1] 

Mandelker N., Dekel A., Ceverino D., Tweed D., Moody C. E., 
Primack J., 2014, MNRAS, 443, 3675 [2.2] 


Martig M., Bournaud F., Teyssier R., Dekel A., 2009, ApJ, 707, 
250 [6.2, 7.1] 

McGrath E. J., Stockton A., Canalizo G., lye M., Maihara T., 
2008, ApJ, 682, 303 [7.1] 

Mitchell P. D., Lacey C. G., Cole S., Baugh C. M., 2014, MNRAS, 
444, 2637 [1] 

Moody C. E., Guo Y., Mandelker N., Ceverino D., Mozena M., 
Koo D. C., Dekel A., Primack J., 2014, MNRAS, 444, 1389 
[ 2 ] 

Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch 
F. C., Maccio A. V., Naab T., Oser L., 2010, ApJ, 710, 903 
[2.2, A] 

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121 
[2.2, A] 

Munoz J. A., Peeples M. S., 2015, MNRAS, 448, 1430 [1] 

Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569 

[7.1] 

Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191 

[ 2 . 2 ] 

Neistein E., Dekel A., 2008, MNRAS, 388, 1792 [1, 4.1, 4.2] 
Neistein E., van den Bosch F. C., Dekel A., 2006, MNRAS, 372, 
933 [1] 

Nelson E., et al., 2014, Nature, 513, 394 [6.1] 

Noeske K. G., et al., 2007a, ApJ, 660, L43 [1, 4.2, 4.3] 

Noeske K. G., et al., 2007b, ApJ, 660, L47 [1] 

Noguchi M., 1999, ApJ, 514, 77 [1] 

Pannella M., et al., 2009, ApJ, 698, L116 [4.2] 

Pannella M., et al., 2015, ApJ, 807, 141 [1] 

Peng Y.-j., et al., 2010, ApJ, 721, 193 [1, 7.1] 

Rees M. J., Ostriker J. P., 1977, MNRAS, 179, 541 [7.1] 

Renzini A., Peng Y.-j., 2015, ApJ, 801, L29 [1] 

Rodighiero G., et al., 2010, A&A, 518, L25 [1] 

Rodighiero G., et al., 2011, ApJ, 739, L40 [1, 4.3] 

Salim S., et al., 2007, ApJS, 173, 267 [1] 

Salim S., Fang J. J., Rich R. M., Faber S. M., Thilker D. A., 2012, 
ApJ, 755, 105 [7.1] 

Sargent M. T., et al., 2014, ApJ, 793, 19 [1, 6.1] 

Schawinski K., et al., 2014, MNRAS, 440, 889 [7.1] 

Schreiber C., et al., 2015, A&A, 575, A74 [1, 3, 4.2, 4.3, 6.2] 
Scoville N., et al., 2015, preprint, (arXiv: 1505.02159) [1, 6.1] 
Silverman J. D., et al., 2015, ApJ, 812, L23 [1, 6.1] 

Snyder G. F., et al., 2015, MNRAS, 454, 1886 [2] 

Sparre M., et al., 2015, MNRAS, 447, 3548 [1] 

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 
2014, ApJS, 214, 15 [1, 3, 4.2, 4.3] 

Stark D. P., Schenker M. A., Ellis R., Robertson B., McLure R., 
Dunlop J., 2013, ApJ, 763, 129 [4.2] 

Szomoru D., Franx M., van Dokkum P. G., 2012, ApJ, 749, 121 

[7.1] 

Tacchella S., Trent! M., Carollo C. M., 2013, ApJ, 768, L37 [1] 
Tacchella S., Dekel A., Carollo C. M., Ceverino D., DeGraf C., 
Lapiner S., Mandelker N., Primack J. R., 2015a, preprint, 
(arXiv:1509.00017) [2, 5.3, 5.4, 7.1, A] 

Tacchella S., et al., 2015b, Science, 348, 314 [5.3, 6.2, 7.1, 7.2] 
Tacchella S., et al., 2015c, ApJ, 802, 101 [7.1] 

Tacconi L. J., et al., 2010, Nature, 463, 781 [1, 6.3.1] 

Tacconi L. J., et al., 2013, ApJ, 768, 74 [1, 3, 6.3.1] 

Thompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167 

[ 2 . 1 ] 

Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V., 
Hernquist L., 2014, MNRAS, 438, 1985 [1] 

Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., 
Dekel A., 2002, ApJ, 568, 52 [1] 

Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012, 
ApJ, 754, L29 [1, 4.3] 

Whitaker K. E., et al., 2014, ApJ, 795, 104 [1, 6.2] 

Williams R. J., Maiolino R., Santini P., Marconi A., Cresci G., 
Mannucci F., Lutz D., 2014, MNRAS, 443, 3780 [6.1] 


MNRAS 000, 1-27 (2016) 


Compaction, Depletion and Replenishment on the Star-Forming Main-Sequence 23 



logloMvir[IVIo] 


Figure Al. Ratio of stellar mass to halo mass as a function of 
halo mass. The red, green and blue point show our simulations 
at 2 = 2, 1.5 and 1.0, respectively. The large symbols show the 
median (at its 1 a scatter) at these three redshifts. The horizontal 
black dashed line shows the cosmic baryon fraction of the Uni¬ 
verse. The solid lines show the (Behroozi et al. 2013) relations. 
Our simulations lie a factor of 2 — 6 above the Behroozi et al. 
(2013) relation at 2 ~ 1 — 2, but are in agreement with recent ob¬ 
servations of SFGs at 0.5 < 2 : < 2.6 (Burkert et al. 2015), which 
are shown as stars and their median as dashed line. 
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APPENDIX A: STELLAR-TO-HALO MASS 
RELATION 

One of the key quantity for simulations to match is the stel¬ 
lar mass made within a given dark-matter halo. In this sec¬ 
tion, we compare the stellar-to-halo mass (M* — Mvir) rela¬ 
tion of the simulations with results from abundance match¬ 
ing (Conroy & Wechsler 2009; Moster et al. 2010, 2013; 
Behroozi et al. 2010, 2013) and observed kinematics (Burk¬ 
ert et al. 2015). We have already presented a similar com¬ 
parison in our companion paper (Tacchella et al. 2015a). 

Specifically, we compare our simulations with data from 
Behroozi et al. (2013) and Burkert et al. (2015). Behroozi 
et al. (2013) derive the M* — Mvir relation by first assuming 
a general form of this relation, and then obtains the param¬ 
eters by simultaneously performing abundance matching of 
the mass functions at different redshifts and comparing the 
consequent information on star formation with a variety of 
observational data, such as the evolution of the MS and the 
global SFR density of the Universe. Burkert et al. (2015) 
derive dark-matter halo masses from recent observations of 
the Ha kinematics of « ~ 0.8 — 2.6 galaxies, i.e., these es¬ 
timates are totally independent of abundance matching. To 


derive the halo mass, they assume implicitly that the spe¬ 
cific angular momentum of the baryons on the scale of the 
dark halo is the same as that of the dark matter component. 
Hence, their results depend on the angular momentum dis¬ 
tribution on the halo scale (of both baryons and dark mat¬ 
ter), as well as on any re-distribution of angular momentum 
between different baryonic components (e.g., inner and outer 
disk, outflow, bulge). 

Figure Al shows the ratio of stellar mass to halo mass as 
a function of halo mass for our simulations and the estimates 
from observational data (Behroozi et al. 2013; Burkert et al. 
2015). At z = 1, the simulated galaxies have a median halo 
mass of logjQ Mvir = 11-7 ± 0.3 and a stellar to halo mass 
ratio of logj^Q M*/Mvir = —1.3 ± 0.2. Burkert et al. (2015) 
found logj^Q M*/Mvir = —1.5 ± 0.3 at logjo Mvir = 11.7, 
which is consistent with our simulations, while Behroozi 
et al. (2013) found logj^Q M*/Mvir = —1.8 ± 0.2 via abun¬ 
dance matching, which is a factor of 3 lower than our simu¬ 
lated estimate. At z = 2, we find only little evolution in the 
halo-to-stellar mass ratio: our simulated galaxies have a me¬ 
dian halo mass of logjg Mvir = 11.5±0.3 and a stellar to halo 
mass ratio of logjg M*/Mvir = —1.4 ± 0.2. Via abundance 
matching, Behroozi et al. (2013) found logj^g M*/Mvir = 
—2.3 ± 0.2 at logiQ Mvir = 11.5, but Burkert et al. (2015) 
found logjQ M*/Mvir = —1.6±0.3, which is again consistent 
with our simulations. 

We conclude that our simulations produce stellar to 
halo mass ratios that are in the ballpark of the values esti¬ 
mated from observations, and within the observational un¬ 
certainties. Therefore, it is sensible and adequate to use our 
simulations with the adopted feedback prescription, while 
bearing in mind the factor-of-two uncertainties. 


APPENDIX B: ANALYSIS OF THIN 
TIMESTEPS 

For six galaxies (11, 12, 14, 25, 26, and 27), there are thinner 
time-steps available for the analysis. Our standard temporal 
resolution is Aa = 0.01), which corresponds to ~ 100 Myr 
between each snapshot. This is roughly half the orbital time 
at the disc edge, and therefore it should be short enough to 
trace galaxy internal processes. 

We nevertheless compare our obtained results with this 
standard temporal resolution with the high resolution. The 
thinner time-steps are separated by Aa = 0.0005 — 0.0007, 
they have a ~ 20-times higher resolution. In Figure Bl, we 
plot the evolution along the MS (distance from the MS ver¬ 
sus total stellar mass) for the higher (upper panels) and stan¬ 
dard (bottom panels) resolution time-steps. We find that the 
standard resolution tracks the higher resolution very well: 
all the main features over long (~ 1 Gyr) as well as short 
(~ 100 Myr) are there. Only the very short term fluctua¬ 
tions (< 100 Myr) are underestimated, leading to a deficit 
of power of short-term changes (see also Section 5.5). 


APPENDIX C: GRADIENTS ACROSS THE MS 

One of the main goals of this paper is to determine the 
galaxy properties across the MS. As discussed in Section 5, 
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Figure Bl. Evolutionary tracks about the universal MS comparing short and standard timesteps between snapshots. We compare for 
six galaxies the higher resolution time-steps (top panels, snapshots separated by Aa ^ 0.0006) with the standard resolution time-steps 
(bottom panels, Aa Ri 0.01). We find that the standard resolution tracks the main features of the evolution along the MS rather well. 


if we want to determine these MS gradients, we have to take 
into account that galaxy properties may also evolve with 
cosmic time. We therefore have to correct for the systemic 
time evolution of galaxies near the MS ridge. 

Cl Gas-related gradients 

Figure Cl is a more extended version of Figure 7. We again 
investigate the following four key quantities are investigated 
from top to bottom: the central gas density within 1 kpc 
(Pgas,i), the total gas mass (Mgas), gas to stellar mass ratio 
(fgs), and depletion time (tdep)- 

The left column shows these (uncorrected) quantities 
as a function of the distance from the MS ridge (Ams) for 
all galaxies at 2 ; = 1 — 6. The colour coding of the points 
corresponds to redshift. We confirm the trends found in Sec¬ 
tion 5.2: Pgas,i increases systematically by about one order 
of magnitude from below to above the main sequence, i.e., as 
Ams varies from —0.5 to -1-0.5 dex. The trend for the total 
Mgas is weaker, partly contaminated by high-redshift galax¬ 
ies with low gas mass near the MS ridge, fgs shows not only 
a correlation with Ams, but also a clear redshift evolution 
which has been highlighted also in Figure 1. The depletion 
time tdep shows a negative correlation with Ams, with sev¬ 
eral outliers with high tdep at 2 > 3. Those points indicate 
the quenching attempts which were not successful, due to 
the high redshift and a too low halo mass (see Section 6.3.3 
below). 

The best fits for the 2 -dependence (f(z) = ^ ( x 

logj^Q(l-|- 2 )) are shown in the middle-left panels of Figure Cl 
with the colour-coding corresponding to stellar mass M*. 
We find a steep redshift evolution for Mgas and fgs with ( — 
— 1.71 ± 0.15 and -1-1.63 ± 0.10, respectively. For Pgas.i and 
tdep, the 2 -evolution is much shallower with ( = -|-0.30±0.21 
and —0.39 ±0.11, respectively (Table 2). 


In a second step, we correct for the stellar mass depen¬ 
dence, which is shown in the middle-right panels of Fig¬ 
ure 7. All four quantities depend on M*, largely due to 
the strong correlation between gas fraction /gas and M* 
(see Section 3 and Figure 1). We fit the M*-dependence 
of the 2 -corrected quantities with the following relation: 
g{M^) = g+yx (logiQ(M*) —10.5). We find 7 = ±0.36±0.04, 
±0.32±0.04, -0.15T0.02, and -0.19T0.02 for pgas.i, Mgas, 
Mgas/M©, and tdep, respectively (Table 2). 

The right panels of Figure Cl show the key quanti¬ 
ties (pgas.i, Mgas, fgs, and tdep) Corrected for their sys¬ 
tematic 2 -evolution f{z) and M*-dependence g{Ms,) (same 
as Figure 7). By construction, we find for all key quanti¬ 
ties a tighter correlation than before the correction. Fitting 
these corrected quantities with Q = a + S x Ams, we find 
S = ±0.83 ± 0.04, 5 = ±0.48 ± 0.03, <5 = ±0.54 ± 0.02, and 
5 = —0.43 ±0.02 for Pgas.i, Mgas, fgs, and tdep, respectively 
(Table 2). Interestingly, the gradient of Pgas.i across the MS 
is the steepest. This demonstrates that galaxies above the 
MS ridge have a significantly higher central gas density than 
galaxies below it: going from Ams = 0.5 dex below the MS 
ridge to 0.5 dex above it, we find that the Pgaa.i increases 
by nearly an order of magnitude. 


C2 Stellar-structure gradients 

Following the same approach as before for the gas-related 
gradients across the MS, we focns here on the stellar- 
structure gradients. Figure C2 is a more extended version 
of Figure 8 , showing the central stellar mass density within 
1 kpc (Em*,i, the half-mass radius (Re), and the S&sic in¬ 
dex (n). From the left to the right, we show the uncorrected 
quantities, the 2 -dependence, the M*-dependence, and the 
corrected quantities. The best-fitting values are given in Ta¬ 
ble 2. The main difference to the gas-related gradients is that 
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Figure Cl. Galaxy properties across the MS. Shown from top to bottom are the gas mass within 1 kpc, total gas mass, gas mass 
to stellar mass ratio, and depletion time. In the left panels, we plot these quantities as a function of the distance from the MS Ams- 
The redshift dependence is colour coded. The middle-left panels show the intrinsic redshift dependencies, colour-coded by stellar mass. 
The green line indicates the best-fitting redshift dependence /(z). The middle-right panels show the redshift-corrected quantities as a 
function of stellar mass. The green line indicates the best-fitting mass dependence g{z). The right panels display the quantities corrected 
for redshift and mass (same as Figure 7). The green line indicates the best-fitting, while the red dashed line shows the best-fit of the 
observations by G15. 


there is no correlation between the stellar-structure quanti¬ 
ties and the distance to the MS, and the inferred gradients 
are much shallower. 


APPENDIX D: INFLOW RATE, OUTFLOW 
RATE AND STAR-FORMATION RATE ALONG 
THE MS 

In Section 5.4 we discuss the driver of the MS gradients, 
focusing on the balance between the gas inflow rate (input 


term) and SFR plus gas outflow rate (drainage term). Fig¬ 
ure 9 shows the relation between the rate of change of dis¬ 
tance from the MS, Ams, and the balance between the input 
term and drainage terms of gas in the central 5 kpc, namely 
logj^Q Rskpc = log[inflow rate/(SFR-|-outflow rate)]. Here, 
we split the balance term log^Q Bskpc into its constituents, 
namely inflow rate, outflow rate, and SFR. 

Figure D1 shows the change of the distance from the 
MS, Ams, as a function of inflow rate, outflow rate, and 
SFR, each measured within the central 5 kpc. We hud that 
the individual rate are correlated to a lesser degree than the 
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Figure C2. Same as Figure Cl, but for the stellar structure quantities central stellar mass density within 1 kpc half-mass radius 

(Rc), and Sersic index (n). 
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Figure Dl. The origin of the evolution across the MS. The three panels show Ams a function of the inflow rate, outflow rate, and 
SFR from left to right, i.e., we have split the balance term log^Q ^skpc of Figure 9. The colour coding corresponds to the position on 
the MS, Ams- that the outflow rate and SFR are moderately correlated with Ams? whereas the inflow rate is not correlated 

with Ams- 
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combined balance term Bskpc- The outflow rate and SFR are 
moderately correlated with Ams (rpearson = —0.32), while 
the inflow rate is not correlated with Ams- 
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